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THE  TRANSFER  HAMILTONIAN:  A TOOL  FOR  LARGE  SCALE,  ACCURATE, 
MOLECULAR  DYNAMICS  SIMULATIONS  USING  QUANTUM  MECHANICAL 

POTENTIALS 

By 

DeCarlos  E.  Taylor 
May  2004 

Chair:  Rodney  J.  Bartlett 
Major  Department:  Chemistry 

The  accuracy  of  a molecular  dynamics  simulation  is  contingent  upon  the  quality  of 
the  potentials  used  to  compute  the  intermolecular  forces  on  the  constituent  atoms.  The 
most  accurate  potentials  are  those  which  result  from  first  principles  quantum  mechanics 
computations  such  as  coupled  cluster  theory;  however,  for  large  scale  molecular 
dynamics  which  may  be  on  the  order  of  1000’s  of  atoms,  this  is  computationally  not 
feasible  as  most  ab  initio  potentials  are  limited  to  10’s  or  100’s  of  atoms.  In  order  to 
avoid  this  computational  bottleneck,  the  concept  of  a transfer  Hamiltonian  is  introduced. 
The  transfer  Hamiltonian  is  a low  rank,  quantum  mechanical  operator  whose  matrix 
elements  are  computed  from  parametric  functional  forms.  The  parameters  which 
constitute  the  transfer  Hamiltonian  are  fit  to  ab  initio  coupled  cluster  theory  results 
(ionization  potentials,  geometry,  dissociation  energies,  etc.)  for  several  small  model 
systems  which  mimic  the  local  structure  of  the  target  simulation  system.  The  transfer 


Hamiltonian  is  applicable  to  systems  far  beyond  the  reach  of  ab  initio  coupled  cluster 
theory,  yet  it  retains  the  accuracy  of  the  first  principles  computation.  Further,  since  the 
transfer  Hamiltonian  can  be  rapidly  evaluated,  it  is  useful  for  on-the-fly  molecular 
dynamics  simulations  and  gives  results  which  are  more  accurate  than  those  obtainable 
using  a classical  potential. 


CHAPTER  1 

OVERVIEW  OF  DISSERTATION 

The  focus  of  this  dissertation  is  the  accurate  description  of  molecules  and 
materials  using  molecular  dynamics  simulations.  Chapter  2 introduces  molecular 
dynamics,  defines  the  term  “predictive  simulation,”  and  details  the  reliance  of  the 
simulation  results  on  the  quality  of  the  potential  used.  Also,  a critique  of  some  existing 
potentials  found  in  the  literature  of  both  classical  and  quantum  forms  is  given  and  the 
necessity  for  improved  potentials  delineated. 

Chapter  3 introduces  the  transfer  Hamiltonian  and  discusses  its  formal  foundation 
in  terms  of  coupled  cluster  and  density  functional  theories.  Following  the  more  formal 
arguments,  the  functional  form  of  the  operator  is  outlined  and  procedures  for  computation 
of  the  matrix  elements  presented.  Representative  examples  of  transfer  Hamiltonians  for 
various  systems  are  given  in  Chapter  4,  and  a presentation  of  the  parameterization 
strategy  is  outlined. 

The  fifth  and  sixth  chapters  contain  parameterizations  of  transfer  Hamiltonians 
applicable  to  silica  and  silica  in  the  presence  of  water  respectively.  Several  tests  used  to 
validate  the  potentials  are  discussed,  and  in  Chapter  5 a prototypical  silicon  oxide 
nanorod  is  fractured  via  molecular  dynamics  using  the  new  quantum  and  standard 
classical  potentials.  Chapter  6 presents  simulations  of  the  nanorod  in  the  presence  of 
water  and  demonstrates  the  experimentally  observable  phenomenon  known  as  hydrolytic 
weakening  whereby  water  causes  a dramatic  decrease  in  the  strength  of  silica. 
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CHAPTER  2 

MOLECULAR  DYNAMICS 

Introduction 

This  chapter  introduces  molecular  dynamics  simulation  and  defines  the  term 
“predictive  simulation.”  Also,  an  overview  of  some  of  the  classical  and  quantum 
potentials  available  in  the  literature,  with  a focus  on  those  applicable  to  silica  and  silicon, 
is  presented.  The  shortcomings  of  the  available  potentials  will  be  discussed  and  the  need 
for  more  chemically  accurate  potentials  will  be  demonstrated. 

Predictive  Simulation 

Molecular  dynamics  simulations  have  contributed  immensely  to  the  understanding 
of  the  chemo-mechanical  behavior  of  large  molecules  and  materials  [1].  Given  a set  of 
initial  conditions  (positions  and  momenta),  and  the  force  of  interaction  among  the 
constituent  particles,  molecular  dynamics  computes  the  phase  space  trajectories  of  the 
particles  which  constitute  the  material  of  interest  [2-4],  Given  these  trajectories,  the 
chemo-mechanical  properties  of  materials  and  their  response  to  external  perturbations  can 
be  studied  at  the  atomistic  level.  The  ability  to  conduct  studies  with  such  coarse  grained 
refinement  allows  researchers  to  elucidate  the  underlying  reaction  mechanisms  which 
determine  macroscopic  properties. 

In  classical  molecular  dynamics,  atomic  positions  are  obtained  by  integrating 
Newton’s  equation  of  motion 


Fk  (t)  = mk  rk  (0  = 


dV  (r) 

drk 


(2.1) 


2 
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where  Fk(t)  is  the  force  acting  on  the  kth  particle  due  to  the  other  N-l  nuclei  at  time  t. 

The  particle  positions  and  mass  are  denoted  by  r and  m respectively  and  the  dots  indicate 
time  derivatives.  The  potential,  V(r),  determines  the  interaction  among  the  particles. 

These  equations  are  typically  solved  using  finite  difference  algorithms  as 
analytical  solutions  are  not  possible  for  all  but  the  most  simple  systems  and  potential 
forms  [5].  The  most  efficient  and  commonly  used  finite  difference  algorithms  used  in 
molecular  dynamics  are  the  Verlet  [6]  and  the  Predictor-Corrector  algorithms  [7].  The 
quality  of  any  simulation  based  on  either  of  these  algorithms  is  intrinsically  based  upon 
the  fidelity  of  the  potential  used  in  the  algorithm.  As  an  illustration  of  this  dependence, 
consider  the  Verlet  algorithm. 

The  Verlet  algorithm  is  based  on  a Taylor  expansion  of  the  position  vector  about 

time  t 

r(t  + AO  = r(t)  + f(A/)  + yr(At2)  + jyr(Af3)  + ...  (2.2) 

Writing  equation  2.2  for  (t  - At)  yields 

r (t  - A t)  = r (t)  - r ( A t)  + -r(Ar) — r‘(  A f 3 ) + ...  (2.2) 

2 3 ! 

Addition  of  equation  2.2  and  2.3  cancels  all  odd  order  terms  and  leaves 

r(t+  A t)  = 2 r (t ) - r(t-  At)  + rAt1  + O (At4)  (2.4) 

Equation  2.4  is  an  example  of  the  working  equation  necessary  to  perform  a molecular 
dynamics  simulation  using  the  Verlet  algorithm.  In  actual  use,  velocities  are  also 
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computed  and  used  in  conjunction  with  equation  2.4,  but  for  our  illustrative  purpose,  the 
form  of  equation  2.4  is  sufficient. 

According  to  equation  2.4,  the  particle  positions  at  time  (t  +At)  can  be  obtained 
from  knowledge  of  the  positions  at  the  current  time  t,  positions  at  time  (t  - At),  and  the 
accelerations  (term  involving  the  second  time  derivative).  The  only  sources  of  error  in 
equation  2.4  and  the  resulting  phase  space  trajectories  which  would  result  from  use  of 
this  particular  algorithm  are  the  intrinsic  instability  of  the  finite  difference  approach,  the 
necessary  truncation  of  the  Taylor  series’  (equations  2.2  and  2.3),  and  the  accelerations 
which  enter  equation  2.4.  The  discretization  error  can  be  controlled  with  the  use  of  small 
time  steps,  At,  in  the  numerical  integration  of  the  equations  of  motion.  In  addition,  the 
truncation  error  of  the  formally  infinite  Taylor  series,  as  shown  in  equation  2.4,  is  fourth 
order  in  At.  Provided  that  the  time  step,  At,  is  chosen  to  be  small,  this  truncation  error 
will  become  infinitesimal.  In  practice,  time  steps  on  the  order  of  femtoseconds  are  used; 
therefore  the  truncation  error  can  be  safely  disregarded  in  the  determination  of  the 
trajectories  via  equation  2.4.  The  appropriateness  of  a chosen  time  step  can  be  tested  by 
running  a simulation  and  comparing  the  results  to  a simulation  where  a smaller  time  step 
is  used.  If  there  is  no  appreciable  change  in  the  results  then  an  appropriate  choice  has 
been  made  for  the  time  step. 

Given  that  the  truncation  error  is  negligible,  the  only  remaining  source  of  error  in 
equation  2.4  is  due  to  the  accelerations.  The  accelerations  are  evaluated  according  to 
equation  2.1  (Newton’s  equation  of  motion).  Examination  of  equation  2.1  shows  that  the 
computed  accelerations  are  directly  contingent  upon  the  potential  used  to  model  the 
interaction  among  the  particles.  Therefore,  in  order  for  the  phase  space  trajectories  which 
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result  from  a simulation  to  be  indisputable,  knowledge  of  an  exact,  or  at  least  chemically 
accurate,  potential  is  required.  Given  such  a potential,  which  must  result  from  a quantum 
mechanical  treatment,  there  would  be  no  formal  source  of  error  (other  than  numerical 
precision)  inherent  in  the  molecular  dynamics  simulation.  A simulation  of  this  quality  is 
considered  to  be  “predictive”. 

Molecular  Dynamics  Potentials 

Classical  Potentials 

One  solution  to  this  problem  is  to  construct  an  empirical  interatomic  potential, 
which  gives  the  energy  E of  a set  of  particles,  as  an  explicit  mathematical  function 
of  the  set  of  particle  coordinates.  If  this  function  is  sufficiently  easy  to  calculate, 
and  if  it  gives  a sufficiently  accurate  description  of  the  real  physical  system  of 
interest,  then  one  can  perform  realistic  calculations  of  the  properties  of  quite  large 
systems. . . .Of  course,  such  an  approach  will  inevitably  involve  a significant  loss 
of  accuracy,  compared  with  ab  initio  calculations.  - [8,  p.  6991] 

Though  the  need  for  ab  initio  quantum  mechanical  potentials  is  recognized,  the 
choice  of  empirical  potential  becomes  a computational  necessity  as  the  evaluation  of  the 
forces  of  interaction  using  an  ab  initio  wavefunction  technique,  for  systems  with  more 
than  a few  hundred  atoms,  is  mostly  impossible  using  modem  computers.  For  systems 
whose  size  is  amenable  to  ab  initio  treatment,  the  computation  is  still  too  time  consuming 
to  be  practical  for  routine  application  in  molecular  dynamics  where  several  thousand  time 
steps,  each  of  which  requires  an  evaluation  of  the  wavefunction,  may  be  necessary  in 
order  to  observe  the  phenomena  of  interest.  Therefore,  in  practice,  most  molecular 
dynamics  simulations  employ  parameterizable  potentials  of  relatively  simple  form  to 
determine  the  variation  of  the  particle  positions  with  time.  Although  these  potentials  can 
be  rapidly  evaluated  and  are  applicable  to  much  larger  scale  systems  than  possible  with 
ab  initio  techniques,  they  are  in  error  in  comparison  to  results  obtained  from  first 
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principles  quantum  mechanics  computations  on  small  clusters;  thus  they  can  not  be 
considered  predictive.  Furthermore,  classical  potentials  can  not  distinguish  between 
radicals  and  ions,  different  spectroscopic  states,  or  be  used  to  properly  describe  bond 
breaking  in  molecules. 

Two  body  potentials 

There  are  many  parameterizable  classical  pair  potentials  available  for  use  in 
simulations,  each  designed  to  suit  the  problem  of  interest.  Once  a choice  of  functional 
form  has  been  made  for  the  potential,  the  parameters  are  fit  to  be  descriptive  of  the 
system  to  which  the  potential  will  be  applied.  As  an  example,  consider  the  classical  pair 
potentials  of  van  Beest  et  al.  (BKS)  [9]  and  Tsuneyuki  et  al.  (TTAM)  [10]  which  were 
designed  to  be  applicable  to  large  polymorphs  of  silica.  The  BKS  and  TTAM  potentials 
are  very  similar  and  have  the  functional  form 


where  the  first  term  models  the  Coulomb  charge-charge  interaction  between  particles  of 
charge  q.  The  remaining  two  terms  are  collectively  called  the  “Buckingham  term”  and 
are  intended  to  model  the  short-range  repulsive  and  dispersion  interactions  [11].  This 
potential  consists  of  1 1 parameters  consisting  of  the  charges  on  silicon  and  oxygen,  and 
the  three  pair  parameters  a,j,  by,  and  Cjj  for  each  pair  interaction  (O-O, Si-Si, Si-O). 

To  determine  the  parameters,  guided  by  the  low  pressure  crystal  structure  of  silica 
which  is  composed  of  comer  sharing  Si04  tetrahedra,  the  BKS  and  TTAM  potentials 
were  fit  to  reproduce  the  ab  initio  Hartree-Fock  potential  energy  surfaces  of  Si04H4  and 
Si04(4-)  + 4e+  clusters  respectively.  The  BKS  parameterization  was  further  augmented 
with  experimental  data  for  the  elastic  constants  and  unit-cell  dimensions  of  quartz.  The 


(2.5) 
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Hartree-Fock  potential  energy  surface  of  the  Si04  cluster  was  that  corresponding  to  Td 
and  D2h  distortions  of  the  oxygen  atoms  from  their  equilibrium  positions.  Examination  of 
figures  in  the  original  TTAM  paper  indicate  that  the  ab  initio  potential  was  sampled  in 
the  region  of  Req  - 0.2  to  Req  + 0.4  Angstroms  and  that  the  chosen  functional  form  of  the 
pair  potential  reproduces  the  shape  of  the  ab  initio  potential  energy  surface  in  that  region. 

Examination  of  the  BKS/TTAM  parameterization  is  a clear  example  of  the 
seemingly  arbitrary  parameter  choice  evident  in  some  pair  potentials.  The  assignment  of 
the  silicon  and  oxygen  charges  (+2.4  and  -1.2  respectively)  by  BKS  and  TTAM  were 
based  on  Mulliken  population  analyses  of  the  Si04  parameterization  cluster  [12].  It  is 
well  known  that  Mulliken  population  analyses  are  very  much  dependent  on  choice  of 
basis  set  and  may  not  be  a reliable  indication  of  atomic  charge  and  a better  choice  would 
have  been  a charge  analysis  based  on  the  use  of  natural  atomic  orbitals  which  are 
contingent  upon  the  electron  density  and  not  the  basis  set,  though  some  effect  will  still  be 
evident  [13-14]. 

Further,  the  ab  initio  surface  was  only  explored  in  regions  close  to  equilibrium 
with  low  level  Hartree-Fock  calculations.  Although  Hartree-Fock  theory  gives 
equilibrium  geometries  which  are  in  good  agreement  with  experiment,  the  BKS  and 
TTAM  potentials  have  been  used  extensively  to  study  fracture  in  bulk  silica  [15-16]. 
Fracture  entails  the  breaking  of  chemical  bonds  and  sampling  of  configurations  far  from 
equilibrium.  Therefore,  a classical  potential,  which  derives  from  only  near  equilibrium 
reference  data,  has  to  be  used  with  caution  when  applied  to  non  equilibrium  phenomena 
such  as  fracture.  The  inability  of  the  BKS/TTAM  potential  form  to  break  chemical  bonds 
is  clearly  demonstrated  in  Figure  2-1  where  the  potential  energy  surface  for  rupture  of  a 
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brigded  Si-0  bond  (see  figure)  in  pyrosilicic  acid  (Si207H6)  at  the  ab  initio  coupled 
cluster  level  of  theory  (CCSD  with  a DZP  basis)  is  compared  to  that  obtained  with  the 
BKS  potential.  The  CCSD  model  has  a demonstrated  ability  to  achieve  chemical 
accuracy  [17]  however;  the  BKS  potential  grossly  underestimates  the  dissociation  energy 
which  results  from  the  quantum  calculation  and  is  also  in  error  in  the  repulsive  region. 

Also  shown  is  a plot  of  the  forces  which  result  from  the  computed  potential 
energy  surface.  These  forces  are  those  which  would  be  used  in  a molecular  dynamics 
simulation  and  the  classical  forces  are  clearly  in  error  in  comparison  to  the  quantum 
mechanical  forces  in  both  the  repulsive  and  bond  break  regions.  Given  that  the  classical 
potential  does  not  reproduce  the  ab  initio  results,  any  molecular  dynamics  simulation 
based  on  these  potentials  cannot  be  considered  predictive. 

Three  body  potentials 

A natural  first  step  in  improving  the  performance  of  pair  potentials  is  to 
evaluate  interactions  which  involve  atomic  triples  in  conjunction  with  effects  from 
atomic  pairs.  Potentials  of  this  type  are  crucial  to  the  proper  description  of  covalently 
bonded  materials  which  present  a difficult  challenge  for  classical  potentials  because 
complex  quantum  effects,  such  as  hybridization  and  charge  transfer,  have  to  be  described 
by  an  effective  interaction  in  which  the  quantum  mechanical  degrees  of  freedom  have,  in 
some  way,  been  integrated  out  [18].  Further,  general  pair  potentials  favor  the  formation 
of  close-packed  structures,  whereas  most  covalently  bound  systems,  such  as  silicon,  have 
more  open  structures  and  require  a proper  description  of  three  particle  interactions  [19]. 

As  an  example,  consider  the  three  body  potential  of  Stillinger  and  Weber  [20]. 
This  potential  has  been  used  extensively  in  the  study  of  bulk  silicon  and  consists  of  a two 
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CCSD/BKS  POTENTIAL  ENERGY  SURFACE  FOR  PYROSILICIC  ACID  DISSOCIATION 


CCSD/BKS  Forces  for  Pyrosilicic  Acid  Dissociation 


R-Re  (Ang) 


—■—CCSD 

—•—BKS 


Figure  2-1.  Comparison  of  ab  initio  coupled  cluster  (CCSD)  potential  energy  surface  and 
associated  forces  versus  those  obtained  with  the  classical  BKS  pair  potential 
for  pyrosilicic  acid  dissociation.  The  reactive  bond  is  indicated  in  the  figure. 


body  contribution,  and  a three  body  term,  V3,  expressed  as  a separable  product  of  radial 


functions  g(r)  and  an  angular  function 


v3(rij  . rik  ) = <?Oy  )g(rik  )(cos  0 jik 


2 


(2.6) 
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where  ©j*  is  the  angle  subtended  by  atoms  j and  k at  the  vertex  i.  The  V3  term  has  a 
minimum  at  0=109.5°;  thus  the  potential  discriminates  in  favor  of  tetrahedral  bond 
angles  which  are  prominent  in  the  most  stable  zero  temperature  and  pressure 
configuration  of  silicon,  the  diamond  lattice  [21]. 

In  addition  to  the  structure  of  the  diamond  lattice,  the  SW  potential  was 
parameterized  such  that  the  melting  point  and  the  liquid  structure  of  the  resulting  melt 
was  in  agreement  with  experiment.  Although  the  equilibrium  properties  of  silicon  were 
well  characterized  with  the  SW  potential,  a study  conducted  by  Dodson  [22]  on  the 
growth  of  silicon  by  the  addition  of  layers  only  a few  atoms  thick  showed  that  the 
potential  was  inadequate  in  regions  away  from  the  ideal  geometry  and  where  the  sp3 
hybridization  of  the  equilibrium  configuration  had  been  destroyed.  Similar  conclusions 
concerning  the  failure  of  the  SW  potential  must  apply  to  fracture  as  well,  where  non  ideal 
local  geometries  will  certainly  occur. 

Conclusions  analogous  to  those  of  Dodson  were  reported  by  Biswas  and  Hamann 
[23]  and  led  to  the  formulation  of  a new  family  of  three  body  potentials  by  Tersoff  [24], 
These  three  body  potentials  are  fundamentally  different  from  those  of  Stillinger  and 
Weber  in  that  the  strength  of  individual  bonds  is  contingent  upon  the  local  environment 
and  they  are  not  implicitly  biased  to  favor  tetrahedral  bond  angles.  Andreoni  and  Pastore 
[25],  also  motivated  by  the  growth  process  of  silicon,  studied  the  applicability  of  the 
Tersoff  potential  to  small  silicon  clusters  that  exhibit  geometries  which  differ  from  the 
local  equilibrium  configurations  of  bulk  silicon.  In  this  work,  the  geometries  and  relative 
energetics  of  SiN  (N=4-10)  clusters  resulting  from  the  Tersoff  potential  were  compared  to 
first  principles  density  functional  theory  calculations  The  results  showed  that  the 
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geometric  topology  of  the  N=4,8-10  clusters  were  predicted  incorrectly  by  the  Tersoff 
potential  and  that  the  energy  ordering  of  different  local  minima  was  not  in  accord  with 
the  ab  initio  potential. 

As  a final  example  of  the  failure  of  the  three  body  potential  form  when  compared 
to  ab  initio  calculations,  consider  the  potential  of  Vashishta  et  al.  [26]  which  was 
parameterized  for  study  of  amorphous  silica.  This  potential  consists  of  a two  body  part 
which  includes  terms  to  model  the  steric  repulsion  of  ions  of  different  size,  long  range 
Coulomb  interaction,  and  a charge-dipole  interaction  to  model  effects  due  to  polarization 
of  anions.  The  three  body  function  is  a separable  product  of  a radial  function  and  an 
angular  term,  scaled  by  a strength  parameter  which  is  contingent  upon  the  angle  (Si-O-Si 
or  O-Si-O)  type.  The  three  body  term,  similar  to  that  of  Stillinger  and  Weber,  is  biased 
towards  angles  which  correspond  to  those  observed  at  equilibrium.  The  parameters  in  the 
potential  were  taken  from  standard  tables  (polarizability  of  oxygen),  empirically  chosen 
based  on  the  authors’  experience,  and  some  were  numerically  fit  to  reproduce  the 
experimental  melting  temperature  and  pressure  of  silica. 

A comparison  of  the  potential  energy  surface  and  associated  forces  for  the  ab 
initio  and  classical  3-body  potential  of  Vashishta  et.  al  is  shown  in  Figure  2-2.  Although 
the  three  body  potential  includes  higher  order  interactions  than  do  the  pair  potentials,  the 
agreement  with  the  ab  initio  results  is  still  poor. 

Quantum  Potentials 

High  level  ab  initio  molecular  orbital  calculations  are  an  ideal  choice  of  potential 
for  use  in  molecular  dynamics  simulations.  Ab  initio  techniques  are  capable  of  properly 
characterizing  the  geometric  and  energetic  topology  of  systems  both  near  and  far  from 
equilibrium  which,  thus  far,  has  not  been  possible  with  classical  potentials.  Further, 
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Figure  2-2.  Comparison  of  ab  initio  coupled  cluster  (CCSD)  potential  energy  surface  and 
associated  forces  versus  those  obtained  with  the  3-body  potential  for 
pyrosilicic  acid  dissociation.  The  reactive  bond  is  shown  in  Figure  2-1. 

quantum  mechanical  potentials  are  able  to  distinguish  among  different  states  which  is  not 

even  a concept  when  dealing  with  classical  potentials.  This  type  of  flexibility  is 

necessary  when  dealing  with  non-equilibrium  phenomena  such  as  fracture  where  a 

variety  of  dissociation  pathways  are  possible  and  need  proper  treatment.  Although 

quantum  mechanical  potentials  have  some  very  attractive  features,  their  applicability  to 
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large  scale  simulation  is  severely  limited  by  the  computational  cost  of  the  methods.  Due 
to  this  large  computational  scaling,  quantum  mechanical  Hamiltonians  used  in 
simulations  have  been  limited  to  density  functional  and  semi-empirical  forms. 

Density  functional  potentials 

The  most  rigorous  quantum  mechanical  potentials  routinely  used  in  molecular 
dynamics  simulations  are  those  based  on  density  functional  theory  (DFT)  [27].  Although 
density  functional  methods  give  accurate  geometries  and  forces,  the  system  size  (on  the 
order  of  100’s  of  atoms)  and  time  scale  for  simulations  (a  few  picoseconds),  is  limited  by 
the  use  of  these  potentials. 

Semiempirical  potentials  - Tight  binding  and  NDDO  models 

Due  to  the  computational  limitations  of  density  functional  potentials,  semi- 
empirical quantum  mechanical  Hamiltonians,  which  encompass  tight  binding  and  the 
neglect  of  diatomic  differential  overlap  (NDDO)  models,  have  been  used  extensively  in 
molecular  dynamics  simulations.  In  these  schemes,  the  computationally  expensive 
integrals  present  in  ab  initio  methods  are  replaced  by  simple  parametric  forms. 

There  are  many  tight  binding  models  available  in  the  literature  and  they  can  be 
viewed  as  reparameterized  Huckel  methods  which  include  short  range  repulsive  terms 
[28],  The  most  commonly  used  and  computationally  simple  forms  assume  an  orthogonal 
basis  where  each  atomic  center  is  associated  with  a minimal  set  of  valence  atomic 
orbitals.  The  tight  binding  Hamiltonian  (HTb)  is  limited  to  one  and  two  center 
interactions  whose  values  are  computed  in  terms  of  simple,  parameterizable  functions 
[29-31], 

In  the  orthogonal  tight  binding  method,  once  the  Hamiltonian  matrix  has  been 
built,  the  one  particle  states  are  determined  in  a non  self-consistent  manner  by  solving 
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Htb^  i=e^i  (2.7) 

and  the  total  energy  is  given  by 

ETOT  = ^ (0i£\  + ^ Vrep  ( rij  ) (2-8) 

occ  i i#y 

where  (Dj  is  the  orbital  occupation  number  and  the  second  term,  Vrep,  models  the 
electrostatic  repulsion  between  atomic  cores  [31],  Classically,  this  term  is  given  as 
ZaZb/Rab;  however,  in  tight  binding  models,  this  term  is  also  treated  as  a parameterizable 
function. 

The  non  self-consistent  nature  of  the  tight  binding  method,  combined  with  the 
simplicity  of  the  Hamiltonian  matrix  elements,  renders  tight  binding  faster  and  applicable 
to  larger  size  systems  than  first  principles  techniques.  In  self-consistent  methods  (i.e., 
Hartree-Fock)  where  the  Hamiltonian  depends  upon  the  electron  density,  repeated 
diagonalizations,  which  scale  cubically  with  the  size  of  the  basis  set,  are  required.  This  is 
often  the  rate-determining  step  in  the  computation.  However,  the  tight  binding 
Hamiltonian  does  not  depend  on  the  electron  density;  therefore  the  one  particle  states  are 
obtained  after  a single  matrix  diagonalization.  Although  the  non  self-consistent  nature  of 
the  tight  binding  Hamiltonian  is  advantageous  in  regard  to  computational  speed,  it  suffers 
a serious  deficiency  in  that  system  charge  is  ignored  [29].  Computationally,  the  tight 
binding  Hamiltonian  predicts  the  same  eigenvalues  and  eigenvectors,  irrespective  of  the 
number  of  electrons  in  the  system.  As  a result  of  this,  the  only  quantity  which 
distinguishes  states  of  a different  charge  is  the  total  energy  which  is  contingent  upon  the 
orbital  occupation  numbers  (see  equation  2.8).  However,  since  tight  binding  predicts  the 
same  orbital  energies  for  all  charge  states,  the  total  energy  (and  associated  forces) 
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computed  via  equation  2.8  must  be  in  error  since  adding  or  removing  particles  changes 
the  Coulomb  interactions  among  the  electrons  and  will  cause  a shift  in  the  eigenvalue 
spectrum. 

Given  the  wide  variety  of  tight  binding  parameterizations  available  in  the 
literature,  it  is  difficult  to  generalize  the  success  of  the  method  when  compared  to  ab 
initio  results.  Panda  et  al.  [30]  did  a comparison  of  DFT  and  ab  initio  Hartree-Fock 
theory  to  the  orthogonal  tight  binding  parameterizations  of  Kwon  et  al.  [32]  and  Lenosky 
et  al.  [33].  In  this  work,  a comparison  of  binding  energies,  equilibrium  geometries  and 
dipole  moments  was  made  for  small  Siw  clusters  up  to  N=19.  For  the  binding  energies,  it 
was  reported  that  the  parameterization  of  Kwon  et  al.  agreed  with  the  Hartree-Fock 
results  for  clusters  up  through  N=10;  however,  the  Lenosky  et  al.  parameterization 
predicted  binding  energies  which  were  too  low.  Geometrically,  the  tight  binding  models 
used  in  this  work  gave  errors  in  bond  lengths  of  =0.2  Angstroms  in  comparison  to  the 
density  functional  geometries  and  in  some  cases  tight  binding  predicted  a different 
geometric  topology  (D3h  versus  C3V)  than  did  the  density  functional  method.  Overall,  the 
agreement  with  ab  initio  results  for  these  models  is  fair  but,  unfortunately,  this  study  and 
work  similar  to  it  done  by  other  researchers  using  different  tight  binding 
parameterizations,  do  not  sample  non-equilibrium  configurations.  Similar  to  the  situation 
with  classical  potentials,  tight  binding  has  been  used  extensively  to  study  fracture  where 
non-ideal  geometries  will  occur.  However,  the  fidelity  of  the  method  at  configurations 
away  from  ideality  has  not  been  well  reported.  There  have  been  some  studies  of  defect 
sites  (removal  of  an  atom)  in  bulk  silicon  [31],  but  these  studies  are  focused  primarily  on 
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energetic  effects  and  not  the  local  geometry  at  the  defect,  which  is  indicative  of  forces, 
the  determining  factor  for  a predictive  simulation. 

In  addition  to  tight  binding  models,  semi-empirical  Hamiltonians  based  on  the 
neglect  of  diatomic  differential  overlap  (NDDO)  approximation  have  been  used  to 
generate  forces  in  molecular  dynamics  simulations  [34],  The  NDDO  based  methods  are, 
similar  to  tight  binding  methods,  parameterizable  Hamiltonians  limited  to  two  center 
interactions  at  most,  where  the  Coulombic  repulsion  between  atomic  cores  is  included  in 
an  ad  hoc  way  (i.e.,  not  ZaZb/Rab)-  However,  the  NDDO  methods  are  a step  towards  ab 
initio  theory  because  non-zero  one  and  two  electron  integrals  are  computed  explicitly  in 
terms  of  parameterizable  functional  forms  [35], 

Computationally,  NDDO  methods  are  single  determinant,  self-consistent,  Hartree- 
Fock  level  calculations,  where  the  elements  of  the  Fock  matrix  are  built  from 
parameterized  integrals.  Although  electron  correlation  may  be  considered  explicitly, 
typically  via  configuration  interaction,  it  is  generally  regarded  that  NDDO  models 
include  correlation  at  the  Hartree-Fock  level  via  inclusion  of  parameterized  integrals 
which  are  fit  to  experimental  or  correlated  ab  initio  results  [36].  Therefore,  since 
correlation  is  included  implicity  in  the  parameterization,  to  avoid  any  possible  “double- 
counting” of  correlation  corrections,  most  NDDO  results  are  reported  at  the  Hartree-Fock 
level.  The  self-consistent  nature  of  the  NDDO  Hamiltonian  renders  it  applicable  to  both 
open  and  closed  shell  systems  which  is  a distinct  advantage  over  the  tight  binding 
models.  However,  because  of  the  requisite  matrix  diagonalizations  necessary  to  achieve 
self-consistency  of  the  density  and  the  computation  of  two  electron  integrals  included  in 
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the  model,  the  NDDO  methods  are  an  order  of  magnitude  slower  than  tight  binding 
formulations  in  regard  to  computational  speed. 

There  are  many  parameterizations  of  the  NDDO  Hamiltonians  available  for 
routine  use  [37],  Although  NDDO  methods  are  used  to  study  many  properties  such  as 
ionization  potentials,  heats  of  formation,  and  dipole  moments,  the  ability  to  predict 
accurate  geometries  is  most  important  for  molecular  dynamics  because  accurate 
geometries  indicate  an  accurate  treatment  of  forces.  For  ground  state  equilibrium 
geometries,  all  methods  perform  well  with  average  errors  in  bond  lengths  and  angles  of 
=0.05  Angstroms  and  4°  respectively  as  compared  to  experiment,  but  there  are  some 
cases  where  the  methods  produce  much  larger  than  average  errors  [35], 

Although  the  NDDO  methods  are  attractive  models  for  use  in  molecular  dynamics 
simulations  they,  similar  to  other  classes  of  potentials,  are  parameterized  to  reproduce 
equilibrium  properties.  Therefore,  the  potential  energy  surface  in  regions  away  from 
equilibrium  cannot  be  expected  to  be  well  characterized  with  these  models.  As  a 
demonstration  of  the  failure  of  the  NDDO  potentials  at  non  equilibrium  configurations, 
consider  Figure  2-3  which  compares  results  from  ab  initio  CCSD  calculations  to  those  of 
two  existing  NDDO  Hamiltonians,  AMI  and  MNDO/d.  There  it  is  demonstrated,  most 
notably  for  the  forces,  that  the  NDDO  Hamiltonians  fail  to  accurately  describe 
pyrosilicic  acid  in  regions  away  from  equilibrium. 

Conclusion 

This  chapter  introduced  the  term  “predictive  simulation”  and  demonstrated  the 
lack  of  accuracy,  as  compared  to  ab  initio  theory,  for  some  of  the  commonly  used 
classical  and  quantum  potentials.  Given  the  demonstrated  necessity  of  accurate,  rapidly 
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Figure  2-3.  Comparison  of  ab  initio  coupled  cluster  (CCSD)  potential  energy  surface  and 
associated  forces  versus  those  obtained  with  the  semi-empirical  AMI  and 
MNDO/d  potentials  for  pyrosilicic  acid  dissociation.  The  reactive  bond  is 
shown  in  Figure  2-1. 

calculable  potentials  for  molecular  dynamics,  the  concept  of  the  transfer  Hamiltonian  is 
introduced  in  the  following  chapter.  It  will  be  seen  that  the  transfer  Hamiltonian  is  a 
parameterizable  quantum  mechanical  operator  designed  to  accurately  reproduce  ab  initio 
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forces,  the  most  important  quantity  for  accurate  simulations,  yet  the  computational 
scaling  is  much  lower  than  first  principles  calculations.  The  transfer  Hamiltonian  will  be 
shown  to  be  a practical  tool  for  molecular  dynamics  simulation. 


CHAPTER  3 

TRANSFER  HAMILTONIAN 

Introduction 

This  chapter  will  present  two  formal  justifications  of  the  transfer  Hamiltonian  in 
terms  of  coupled  cluster  theory  and  from  a density  functional  viewpoint.  It  will  be  seen 
that  the  transfer  Hamiltonian  can  be  viewed  as  approximations  to  either  of  these  models 
which,  in  principle,  provide  ab  initio  analogs  to  the  parameters  which  comprise  the  semi- 
empirical  transfer  Hamiltonian.  Following  the  more  formal  arguments,  the  functional 
form  of  the  transfer  Hamiltonian  will  be  presented  in  detail. 

Theoretical  Formulation 

It  has  been  shown  that  both  classical  and  semi-empirical  potential  forms  (in  their 
standard  formulations)  are  inadequate  when  non-equilibrium  areas  of  the  potential  energy 
surface  are  sampled.  In  simulations  of  fracture,  where  bonds  are  being  broken,  it  is 
necessary  to  have  a quantum  mechanical  potential  which  is  capable  of  accurately  treating 
both  equilibrium  and  non  ideal  configurations  if  the  simulation  is  to  be  considered 
predictive.  However,  the  potential  must  be  rapidly  calculable  and  applicable  to  systems 
large  enough  to  properly  simulate  the  phenomena  of  interest. 

As  a means  to  this  end,  the  concept  of  the  transfer  Hamiltonian  is  introduced.  The 
transfer  Hamiltonian  is  a low  rank,  self-consistent,  quantum  mechanical  operator  whose 
matrix  elements  are  given  in  terms  of  parameterizable  functions.  The  functional  form  of 
the  transfer  Hamiltonian  may  be  chosen  freely;  however  the  existing  NDDO 
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Hamiltonians  are  a judicious  choice  due  to  their  demonstrated  ability  to  reproduce  ab 
initio  geometries  in  equilibrium  regions.  In  general,  once  a choice  of  functional  form  has 
been  made,  the  parameters  of  the  transfer  Hamiltonian  are  fit  using  numerical 
optimization  algorithms  such  that  the  forces  and  other  quantities  of  interest  given  by  the 
transfer  Hamiltonian  reproduce  those  given  by  ab  initio  theory. 

The  transfer  Hamiltonian,  constructed  in  this  way,  is  a useful  tool  for  molecular 
dynamics  simulations  because  it  provides  a convenient  route  for  the  predictive  simulation 
of  materials  without  the  computational  bottleneck  of  ab  initio  potentials.  Further,  the 
transfer  Hamiltonian  is  superior  to  classical  potentials  and  tight  binding  models  due  to  its 
demonstrated  ability  to  properly  describe  different  charge  states  and  a wide  range  of 
geometric  topologies. 

The  functional  form  of  the  transfer  Hamiltonian  may  be  chosen  freely,  however, 
practitioners  of  semi-empirical  quantum  chemistry,  in  its  NDDO  formulation,  have 
developed  a variety  of  approximate  Hamiltonians  which  are  used  here  as  models  though 
some  variations  are  made  from  the  standard  representations.  Operationally,  the  transfer 
Hamiltonian  calculation  is  simply  an  independent  particle,  self-consistent  field 
computation  based  on  a single  determinant  reference.  However,  due  to  the 
parameterization,  the  energies  and  forces  are  formally  exact.  This  type  of  structure  is 
similar  to  that  of  density  functional  theory  where  the  energies  and  associated  forces  are 
exact  although  the  effective  Hamiltonian  used  in  the  model  is  of  one-particle  type  [38]. 
The  formal  analysis  of  the  transfer  Hamiltonian  can  be  presented  in  terms  of  coupled 
cluster  or  density  functional  theory,  both  of  which,  in  principle,  are  ab  initio 
representations  of  the  transfer  Hamiltonian. 
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Coupled  Cluster  Theory  Formulation 


In  coupled  cluster  theory,  the  ground  state  wavefunction  is  given  by 


V *)=  ‘T  \fo) 


(3.1) 


where  (J)0  is  a single  Slater  determinant  and  T is  an  excitation  operator  which  creates 
single,  double,  up  to  n-tuply  excited  configurations  from  the  reference  determinant.  T 
can  be  written  as 


where  the  standard  creation/annihilation  operators  have  been  used  [39],  A coupled 
cluster  calculation  reduces  to  a determination  of  the  t amplitudes  present  in  each  level  of 
excitation  and,  as  equation  3.2  shows,  it  is  not  based  on  an  independent  particle  treatment 
of  the  electrons  as  is  the  transfer  Hamiltonian. 

The  connection  of  the  independent  particle  transfer  Hamiltonian  to  many  body 
coupled  cluster  theory  can  be  made  via  consideration  of  the  effective  one-particle 
Hamiltonian  which  arises  in  Brueckner-Coupled  Cluster  (BCC)  theory  [40],  The 
Brueckner  model  is  a version  of  coupled  cluster  theory  where  the  orbitals  are  optimized 
so  that  the  single  excitation  amplitudes  are  equal  to  zero  (Ti=0  in  equation  3.2)  In  BCC, 
the  total  energy  is  the  sum  of  the  independent  particle  reference  determinant  energy  and  a 
correlation  term: 


E = {®\e-THeT\<5>)=Eref+Ei 


corr 


(3.3) 


(3.4) 


(3.5) 
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where  the  equations  have  been  written  in  the  spin  orbital  basis  and  i,j  (a,b)  refer  to 
occupied  (virtual)  spin  orbitals.  In  the  following,  indices  p,q  denote  orbitals  in  both 
occupied  and  virtual  spaces,  h is  the  one  electron  Hamiltonian  which  contains  the  kinetic 
energy  and  nuclear-electron  attraction  terms  and  the  fock  matrix  elements  f which  appear 
in  equation  3.4  are  given  by 

fpq  = hpq+^J  ( Pj  ||  Vi)  (3-6) 

i 

The  total  energy  expression  in  equation  3.3  can  be  written  in  terms  of  an  effective 
Hamiltonian,  G: 


E = ~'L(K  + G,i)  (3.7) 

where  the  general  matrix  element  of  the  occupied-occupied  block  of  G is  given  by 

Gik  = fik  + (ab  || kj)ta..  (3.8) 

jab  V 

and  the  occupied- virtual  block  of  the  effective  matrix  is  chosen  to  be  the  Ti  amplitude 
equation 


71  =G;a  =/,«, + (3.9) 

jb  jbc  jkb 


Solution  of  the  set  of  one-particle  like  equations 


G </>  , = e , ,</>  , (3.10) 

yields  a set  of  orbitals  (Brueckner  orbitals)  which  satisfy  Ti=Gia=0  i.e  the  Brueckner 
condition.  These  orbitals  are  then  used  to  generate  a new  set  of  T amplitudes  which 
define  a new  effective  operator  and  the  procedure  is  repeated  until  self-consistency  is 
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achieved.  The  effective  matrix  G and  the  T amplitude  equation  are  solved  repeatedly 
until  the  T amplitudes,  and  consequently  the  G matrix,  do  not  change  in  successive 
iterations. 

The  total  energy  expressed  in  terms  of  the  effective  one-particle  operator  G (see 
3.7),  is  equivalent  in  form  to  that  of  the  independent  particle  Hartree-Fock  model 
(compare  equations  3.4  and  3.7)  and  it  can  be  concluded  that  BCC  is  an  effective  one- 
particle  theory  whose  energy  and  associated  forces  are  exact. 

When  viewed  in  this  way,  the  transfer  Hamiltonian,  Th,  can  be  considered  as  an 
approximation  to  G since  both  operators  yield  energies  and  forces  which  are  formally 
exact.  In  the  Brueckner-coupled  cluster  case,  the  correlation  contribution  to  G is 
included  via  an  iterative  determination  of  the  T2  amplitudes  which  determine  the  matrix 
elements  of  the  effective  one-particle  operator  G.  However,  in  the  Th  case,  the  many 
body  contribution  to  the  energies/forces  is  obtained  via  parameterization  of  the  matrix 
elements.  Since  each  operator  involves  correlation  corrected  matrix  elements  which  yield 
the  exact  energy  and  forces,  G serves  as  an  ah  initio  analog  of  the  semi-empirical 
operator  Th. 

Similar  approaches  have  been  considered  by  Karl  Freed  [41]  and  coworkers  who 
attempt  to  compute  semi-empirical  PPP  Hamiltonian  parameters  from  first  principles.  In 
their  approach,  they  consider  the  matrix  equation 


" PHP 

PH  Q 

-pxp- 

= E 

pw- 

QHP 

QHQ 

QW 

QW 

where  P is  the  primary  space  and  Q is  the  orthogonal  complement  to  P.  Simple  algebraic 
manipulations  starting  from  equation  3.1 1 yield  an  effective  Hamiltonian 


H v = P H P + P H Q (E  - Q H Q )-'  Q H P 


(3.12) 
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where  Hv  acts  only  in  the  primary  space  but  yields  the  exact  energy  E [42-43].  If  the 
primary  space  is  considered  to  be  a ground  state  single  determinant  reference  and  the 
orthogonal  complement  Q includes  all  single,  double,  to  n-tuply  excited  determinants, 
then,  similar  to  the  coupled  cluster  case  (G),  Hv  consists  of  correlation  corrected  integrals 
which  are  ab  initio  analogs  of  semi-empirical  parameters. 

The  coupled  cluster  viewpoint  is  a more  convenient  and  transparent  route  to  an  ab 
initio  representation  of  the  semi-empirical  Hamiltonians  since  it  bypasses  the  energy 
dependence  and  large  matrix  inversion  required  in  the  Hv  approach.  In  practice,  the 
inverse  matrix  is  replaced  by  a perturbation  expansion  and  the  equation  solved  iteratively 
[41].  On  the  contrary,  the  representation  of  the  matrix  elements  of  G are  given  by  closed 
formulas  in  terms  of  many  body  integrals  and  T amplitudes  determined  in  the  coupled 
cluster  calculation  [40].  Further,  for  a given  choice  of  Q space  determinants,  coupled 
cluster  theory  will  include  contributions  due  to  higher  excitation  manifolds  than  the  Hv 
approach  due  to  the  exponential  form  of  the  wave  operator  used  in  coupled  cluster  theory 
[39];  therefore  the  correlation  corrected  integrals  of  G will  contain  more  of  the  inter- 
electron physics  and  provide  a better  estimation  of  the  semi-empirical  parameters  than 
those  of  Hv. 

Density  Functional  Theory  Formulation 

Density  functional  theory  is  an  exact  one-particle  theory  where  it  is  possible,  in 
principle,  to  solve  for  the  exact  energy  (including  electron  correlation)  through  the 
treatment  of  a set  of  one  particle  equations  [38],  The  density  functional  energy  can  be 


written  as 
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^ DFr  l ^ ] E K E E N e \ P \ + E } [ P ] + E x c \ p ] (3.13) 

where  the  final  term,  the  exchange-correlation  energy,  includes  the  stabilizing  effects  of 
electron  exchange  (K  operator  in  Hartree-Fock  theory)  and  correlation  into  a single 
functional  of  the  density.  Given  an  exact  exchange-correlation  energy  functional,  Edft 
would  be  the  exact  energy  although  the  method  is  of  one-particle  type,  analogous  to  the 
transfer  Hamiltonian. 

The  functional  derivative  of  the  exchange-correlation  energy  functional  yields  the 
exchange-correlation  potential.  This  potential,  when  combined  with  the  usual 
expressions  for  the  one  electron  core  Hamiltonian  and  two  electron  coulomb  integrals  (J) 
yields  an  effective  one-particle  operator  G such  that 

G <t)  , = E ,<l>  , (3.14) 

The  resulting  one-particle  functions  {(j)}  yield  a density,  which  when  evaluated  in  equation 
3.13,  render  an  energy  and  forces  which  are  formally  exact.  The  exchange-correlation 
functional,  similar  to  the  role  of  the  parameters  contained  in  the  transfer  Hamiltonian, 
yields  an  effective  operator  that  is  of  one-particle  type  yet  yields  exact  energies  and 
forces,  therefore,  the  transfer  Hamiltonian  can  be  interpreted  as  modeling  the  Kohn-Sham 
operator. 

Strictly  speaking,  in  order  to  proceed  with  this  interpretation,  the  transfer 
Hamiltonian  matrix  elements  should  be  computed  using  the  exchange -correlation 
potential  (or  some  approximation  to  it)  used  to  solve  the  first  principles  Kohn-Sham 
equations.  However,  the  DFT  based  transfer  Hamiltonian  is  treated  via  a HF-DFT 
(Hartree-Fock  DFT)  computation  where  the  general  exchange-correlation  functionals  are 
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applied  to  the  converged  transfer  Hamiltonian  density  determined  at  the  Hartree-Fock 
level  [44],  When  treated  in  this  way,  the  parameters  which  comprise  the  transfer 
Hamiltonian  cannot  be  interpreted  as  modeling  correlation  corrections  as  in  the  coupled 
cluster  case  because  the  correlation  contribution  to  the  energy  is  taken  care  of  by  the 
exchange-correlation  functional.  The  parameters,  in  the  context  of  a transfer  Hamiltonian 
based  on  HF-DFT,  are  simply  approximations  to  the  ab  initio  one  and  two-electron 
integrals. 

In  summary,  the  transfer  Hamiltonian  is  an  effective  one  particle  operator  which 
yields  exact  energies  and  forces  by  virtue  of  parameters.  This  operator  is  similar  to  the 
effective  one-particle  operator  which  arises  in  Brueckner  coupled  cluster  theory  or  the 
Kohn-Sham  operator  of  density  functional  theory,  both  of  which  provide  ab  initio 
representations  of  the  transfer  Hamiltonian. 

Functional  Form 

Given  a coupled  cluster  or  DFT  based  transfer  Hamiltonian,  the  parameterizable 
one  and  two  electron  integrals  over  atomic  basis  functions  have  to  be  computed  in  a 
suitable  fashion  and  the  Hartree-Fock  self  consistent  field  equations  solved.  In  ab  initio 
models,  these  integrals  are  almost  universally  computed  using  gaussian  type  basis 
functions  for  which  the  requisite  poly-centric  integrals  can  be  computed  analytically  [45]. 
However,  for  the  transfer  Hamiltonian,  these  integrals  are  determined  in  terms  of 
parameterizable  functions. 

The  chosen  integral  approximations  in  the  coupled  cluster  or  density  functional 
variant  of  the  transfer  Hamiltonian,  are  taken  from  the  established  field  of  semi-empirical 
theory  in  its  modified  neglect  of  differential  overlap  (MNDO)  formulation  [46]. 
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In  MNDO,  the  zero  differential  overlap  approximation, 

♦>/(!)]  = ^,,/[^,(l),0,(l)]  (3.15) 

is  invoked  where  O^and  Ov  are  basis  functions  on  centers  A and  B respectively.  The 
implication  of  equation  3.15  is  that  all  integrals  which  involve  multi-center  charge 
distributions  are  zero  which  limits  the  Hamiltonian  to  two-center  interactions  at  most. 

MNDO  assigns  to  each  atom  a set  of  valence.  Slater  atomic  orbital  basis  functions 
[47].  Therefore,  each  first/second  row  atom  receives  a set  of  sp  Slater  basis  functions 
(for  a total  of  4)  and  all  higher  row  atoms  receive  a set  of  spd  (9)  orbitals.  In  some  cases, 
a set  of  polarization  functions  may  be  used  for  hydrogen  and  the  first/second  row  atoms. 

In  that  case,  hydrogen  is  assigned  an  additional  set  of  p orbitals  and  the  first/second  row 
atoms  are  associated  with  a set  of  spherical  harmonic  d functions. 

MNDO  is  a valence  electron  model,  therefore  the  effective  nuclear  charge  for  each 
atom  is  given  by  Z’=Z-n  where  n is  the  number  of  core  electrons  for  the  atom.  Given  the 
basis  set  and  nuclear  charges,  the  one  and  two  electron  integrals  are  computed 
parametrically  as  follows  [35]: 

(in  the  following  equations  (p,v)  and  (k,A.)  are  basis  functions  assigned  to  centers  A and  B 
respectively) 

One  Center  - One  Electron  Integrals 

For  basis  functions  p and  v on  atomic  center  A,  the  one  electron  integrals  are 
given  by 

H = U nn  - Z Z B ( M H I s g s B ) (3.16) 

A * B 


for  the  diagonal  elements  and 
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H „ „ = - £ !■'»*»)  (3.17) 

A * B 

for  the  off  diagonal  contributions.  in  equation  3.16  is  treated  as  a parameter  and  is 
interpreted  as  the  one  center  atomic  orbital  energy  of  the  p’th  orbital  (p=s,p,d  orbital). 
This  parameter  can  be  taken  from  atomic  spectral  data  or  varied  in  the  parameterization 
of  the  transfer  Hamiltonian.  If  the  latter  is  chosen,  care  must  be  taken  to  ensure  that  the 
U parameters  for  each  atom  have  the  proper  energy  ordering.  In  other  words,  the  U 
parameter  for  an  s type  basis  function  should  be  more  negative  than  that  for  p and  so  on. 
This  is  consistent  with  the  interpretation  that  U models  the  one  center  atomic  orbital 
energy  and  the  work  required  to  remove  an  electron  from  an  s orbital  is  greater  than  that 
for  a p orbital  and  similarly  for  the  d orbital. 

As  equations  3.16  and  3.17  show,  the  attraction  of  an  electron  with  charge 
distribution  centered  on  atom  A for  a nucleus  B is  modeled  in  terms  of  the  effective 
nuclear  charge  of  center  B,  weighted  by  the  two-electron  repulsion  integral  involving  the 
charge  distribution  of  center  A and  the  spherically  symmetric  distribution  |ss)  of  center  B. 
This  formula  was  determined  empirically,  without  theoretical  justification  [48]. 

Two  Center  - One  Electron  Integrals 

Although  the  ZDO  approximation  implies  that  this  class  of  integrals  is  zero,  they 
are  retained  in  the  model  because  it  was  found,  by  the  early  practitioners  of  semi- 
empirical  methods,  that  bonding  between  atoms  was  much  improved  by  their  presence  in 
the  model  [36],  For  a two  center  charge  distribution  p on  A and  k on  B,  the  two  center 
core  Hamiltonian  matrix  element  is  computed  as 

H H K ~ — M + P K UK 


(3.18) 
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where  S^K  is  the  overlap  integral  between  the  two  basis  functions.  P is  a freely  variable 
atomic  parameter  with  a separate  value  for  each  basis  function  type  (s,p,d)  assigned  to 
each  center.  The  value  of  this  matrix  element  is  also  contingent  upon  the  Slater  orbital 
exponent  which  determines  the  value  of  the  overlap  integral.  The  orbital  exponents  are 
also  considered  to  be  variable  parameters. 

One  Center  - Two  Electron  Repulsion  Integrals 

The  one  center,  two  electron  repulsion  integrals  are  treated  differently  depending 
on  the  type  of  basis  functions  assigned  to  each  center  [48],  For  an  atom  with  an  sp  basis 
set,  there  are  only  6 two  electron  repulsion  integrals  which  are  non-zero  by  symmetry. 

The  six  integrals,  are  given  by: 

G ss  = (s  j | i s ) G sp  = ( s s | /?  p ) 

G pp  = ( P P I P P ) H sP  = ( s p \ s p ) (3.19) 

Gp2  = (pp\p'p')  Hpp  = (pp'\pp') 

For  an  sp  basis  set,  due  to  the  small  number  of  integrals  and  to  add  more  flexibility  to  the 
parameterization,  these  integrals  are  treated  as  variable  parameters.  On  the  contrary,  for  a 
basis  set  which  includes  d orbitals,  due  to  the  large  number  of  non-zero  integrals,  all  the 
one  center  two  electron  repulsion  integrals  are  computed  analytically  in  terms  of  Slater- 
Condon  F/G  factors  [49],  These  integrals  are  computed  with  an  additional  set  of  orbital 
exponents  which  differ  from  the  orbital  exponents  used  to  compute  the  overlap  integrals 
that  determine  the  two-center  core  Hamiltonian  matrix  elements. 

Two-Center  Two  Electron  Repulsion  Integrals 

Due  to  the  ZDO  approximation,  all  three  and  four  center  two  electron  repulsion 
integrals  are  neglected.  The  two  center  integrals  which  remain  are  treated  by  the  point 
charge  model  of  Thiel  and  Dewar  [50].  In  this  model,  the  one  center  charge 
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distributions  on  each  center  are  expanded  in  terms  of  multipoles  and  each  multipole  is 
then  represented  by  an  appropriate  point  charge  configuration  centered  on  each  nucleus. 
Finally,  the  interaction  of  each  point  charge  of  center  A with  those  of  B is  evaluated 
using  the  Ohno-Klopman  formula  and  summed  to  give  the  two  electron  repulsion  integral 
[51]. 

The  value  of  this  class  of  integrals  is  parameterized  in  terms  of  the  orbital 
exponents  assigned  to  each  nucleus.  The  orbital  exponents  determine  the  distance  of 
each  point  charge  from  its  respective  nucleus.  Therefore,  diffuse  orbital  exponents  render 
point  charge  configurations  which  are  greater  distances  from  each  nucleus  and  yields 
larger  values  of  the  two  center  repulsion  integrals.  This  is  similar  to  ab  initio  models 
where  diffuse  orbitals  have  a larger  radial  extent  and  yield  larger  integral  values  [52]. 

In  addition  to  the  orbital  exponents,  the  two-center  repulsion  integrals  also  depend 
on  the  value  of  the  one  center  repulsion  integrals  which  may  also  be  parameters  as 
described  previously.  The  Ohno-Klopman  formula  which  is  used  to  evaluate  the 
interaction  between  point  charges  i and  j is  given  by 

e2[R-j  + ( Po  + Po)2V  (3.20) 

where  R,j  is  the  separation  between  the  point  charges.  The  additive  terms  pA  and  pB  are 
constants  assigned  to  each  atom  which  are  determined  such  that  the  value  of  the 
corresponding  one  center  two  electron  repulsion  integral  is  reproduced  in  the  limit  Rij=0. 
Computationally,  once  the  one-center  integrals  are  determined,  the  equations  which 
determine  the  additive  terms  pA  and  pB  are  solved  numerically  for  each  atom  as  an 
analytic  solution  is  only  possible  for  the  one  center  (ss|ss)  repulsion  integral  [50]. 
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Given  the  point  charge  separations  and  the  additive  terms  which  enter  the  Ohno- 
Klopman  formula  for  each  nucleus,  the  two  electron  integrals  can  be  rapidly  evaluated. 
Due  to  the  large  number  of  these  integrals  (i.e.  2025  integrals  for  two  atoms  with  spd 
basis  sets),  they  are  computed  in  a local  coordinate  system  where  symmetry  reduces  the 
number  of  non-zero  integrals  which  need  be  evaluated.  The  resulting  symmetry  integrals 
in  the  local  frame  are  subsequently  transformed  into  integrals  over  molecular  coordinates 
using  p and  d orbital  rotation  matrices  [53]. 

Nuclear  Repulsion  Energy 

In  ab  initio  theory,  the  nuclear  repulsion  between  two  atomic  cores  is  given  by 
ZaZb/Rab>  however,  in  the  MNDO  model  this  term  is  represented  by  a parameterizable 
function.  The  core-core  repulsion  energy  is  given  by 


where  the  a exponents  of  the  lead  term,  and  the  (a,b,c)  parameters  of  the  second  term  are 
variable  parameters. 

Computational  Procedure 

The  transfer  Hamiltonian  is  merely  a parameterized  Hartree-Fock  computation 
where  the  parameterization  is  in  terms  of  the  integrals  over  the  basis  functions  and  the 
intemuclear  repulsion  energy  expression,  as  presented  in  the  previous  sections.  A 
summary  of  the  available  parameters  is  given  in  Table  3-1.  The  parameterized  integrals 
are  summed  into  matrix  elements  using  the  usual  equations  from  Hartree-Fock  theory  and 
the  total  energy  and  density  are  determined  self-consistently.  In  the  density  functional 
based  transfer  Hamiltonian,  the  converged  density  is  used  to  evaluate  the  chosen 


E(A, B)  = ZAZB (sAsA  | sBsB )(1  + <?  aARAB+e  ccbRab) 


3.21 
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exchange-correlation  functional,  otherwise,  the  electronic  energy  is  computed  as  in 
Hartree-Fock  theory  via  ,5*Tr(P[H+F]) . P is  the  one-particle  density  matrix  and  H/F 
have  their  usual  meanings. 

Table  3-1:  Listing  of  the  atomic  parameters  which  define  the  transfer  Hamiltonian 


Parameter 

Description 

Uss,Upp,Udd 

One  center  atomic  orbital  energies  (eqn.  3.16) 

Zs,Zp,Zd 

Orbital  exponent  for  overlap  integrals  and  charge 
separations 

Zsn,Zpn,Zdn 

Orbital  exponent  for  one  center,  two  electron 
repulsion  integrals  (d  orbital  basis) 

Ps,Pp,pd 

Beta  resonance  parameter  for  core  Hamiltonian 
matrix  elements  (eqn.  3.18) 

Gss,Gsp,Gpp,Gp2,Hsp,Hpp, 

One  center  two  electron  repulsion  integrals  - sp 
basis  (eqn.  3.19) 

a,a,b,c 

Core  repulsion  energy  parameters  (eqn.  3.21) 

Parameterization 

The  parameterization  of  the  transfer  Hamiltonian  is  of  extreme  importance  as  the 
accuracy  of  the  model  is  contingent  upon  the  parameters  which  are  used  to  build  the 
operator.  In  order  for  the  parameterization  to  proceed,  it  is  necessary  to  have  a set  of  ab 
initio  reference  data  which  the  transfer  Hamiltonian  should  reproduce  and  an  algorithm  to 
determine  the  optimal  set  of  parameters.  The  reference  data  can  be  any  quantity  which 
should  be  reproducible  by  the  transfer  Hamiltonian;  however,  in  the  context  of  molecular 
dynamics,  the  forces  are  the  most  important  properties.  The  parameterization  is  certainly 
not  limited  to  forces,  as  one  can  also  fit  ionization  potentials,  dipole  moments,  vibrational 
frequencies,  activation  barriers  or  the  density  itself.  It  should  be  clear  that  the  larger  the 
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target  reference  data  set,  the  worse  the  transfer  Hamiltonian  will  perform,  numerically, 
for  each  property.  This  is  more  a consequence  of  the  convergence  of  the  chosen 
numerical  fitting  algorithm  than  an  inherent  limitation  of  the  Hamiltonian.  Therefore,  it 
is  recommended  to  focus  on  one  or  two  types  of  reference  quantities  (i.e.  energy  and 
forces)  for  each  parameterization.  Given  the  reference  data,  the  parameters  are 
determined  using  a combination  of  genetic  and  Newton-Rapshon  optimization 
algorithms. 

Conclusion 

This  chapter  presented  the  formal  justification  and  the  functional  form  of  the 
transfer  Hamiltonian.  It  was  seen  that  the  integral  approximations  are  simple  in  form  and 
through  their  parameterization,  provide  a convenient  route  for  the  introduction  of  ab 
initio  forces  into  molecular  dynamics  simulations,  provided  the  parameters  can  be 
determined.  That  task  is  discussed  in  the  following  chapter. 


CHAPTER  4 

EXAMPLE  TRANSFER  HAMILTONIANS 

Introduction 

This  chapter  will  detail  the  parameterization  of  the  transfer  Hamiltonian.  The 
performance  of  the  optimization  algorithms  for  this  task  will  be  discussed  and  converged 
parameter  sets  for  two  representative  molecules,  ethane  and  disilane,  will  be  presented  as 
a demonstration  of  the  procedure  and  the  accuracy  of  the  model. 

Parameterization  Procedure 

The  parameters  of  the  transfer  Hamiltonian  are  determined  using  genetic  and  quasi- 
Newton  optimization  algorithms.  The  complexity  of  the  transfer  Hamiltonian  lies  not  in 
the  concept  of  the  method  but  in  the  difficulty  of  the  parameterization.  It  is  often  stated 
that  parameterizing  is  an  art,  not  a science.  The  transfer  Hamiltonian  certainly 
corroborates  that  statement. 

Genetic  Algorithms 

Genetic  algorithms  (GA)  are  stochastic  optimization  methods  which  take  their 
terminology  and  conceptual  basis  from  evolutionary  biology  [54].  GAs  are  based  on 
Darwin’s  theory  of  evolution  where  the  strongest,  most  fit  individuals  survive  and 
produce  offspring  which  have  the  favorable  characteristics  of  the  parent  generation.  The 
idea  is  to  have  a population  of  individuals,  each  characterized  by  a set  of  genes.  The 
most  fit  parent  individuals  are  allowed  to  produce  offspring  by  sexual  reproduction,  and 
new  and  improved  generations  are  produced  while  those  less  fit  individuals  are  removed 
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from  the  population.  After  many  generations,  a population  of  strong  individuals  will  be 
created. 

GAs  are  a computational  manifestation  of  evolutionary  biology.  In  the  context  of 
the  transfer  Hamiltonian,  the  population  is  composed  of  trial  parameter  sets  p,  which  are 
used  to  compute  the  one  and  two  electron  repulsion  integrals  described  in  chapter  3.  The 
fitness  of  each  parameter  set  is  determined  by  computing  the  error,  f , between  properties 
computed  with  the  transfer  Hamiltonian  and  ab  initio  reference  data  using  the  equation, 

/ ( P ) = X I O)  l h - Q)  , R e f | (4.1) 

I 

where  f(p)  implies  the  implicit  dependence  on  the  parameter  set  p and  co™  and  coRef  are 
the  values  of  the  particular  property  computed  using  the  transfer  Hamiltonian  and 
contained  in  the  reference  data  set  respectively.  The  most  fit  (those  with  the  smallest 
error)  parameter  sets  are  allowed  to  produce  new  parameters  in  a way  which  models 
evolutionary  biology,  via  the  exchange  of  genetic  information  between  two  parent 
individuals  through  sexual  reproduction.  Computationally,  sexual  reproduction  is 
modeled  by  treating  the  digits  which  represent  each  parameter  as  genes.  A point  on  each 
gene  is  randomly  selected,  and  all  the  digits  to  the  right  of  the  randomly  selected  point 
are  interchanged,  producing  a new  parameter.  This  procedure  of  fitness  assessment  and 
parameter  reproduction  is  repeated  until  convergence  is  achieved. 

Newton-Raphson  Algorithms 

These  well  known  methods  of  optimization  [55-56]  are  gradient  following  schemes 
where  the  function  to  be  minimized  is  expanded  to  second  order  as: 

f(p)  = f{p0  ) + gT(p-Po)  + -^(P~Po)T  h(P-Po)  (4-2) 
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The  derivative  of  the  error  function,  equation  4.1,  with  respect  to  the  parameters  cannot 
be  done  analytically;  therefore,  numerical  differentiation  with  respect  to  each  parameter 
is  necessary.  Requiring  the  gradient,  g of  equation  4.2  to  be  zero  yields 

(p-Po)  = H-'s  (4.3) 

which  is  a prescription  for  the  generation  of  the  new  parameter  set  p from  the  old  set  p0, 
provided  the  gradient  and  Hessian  are  available.  Algorithms  which  avoid  the  direct 
determination  of  the  Hessian  matrix  by  successive  approximations  are  also  available  [57]. 
The  optimization  proceeds  with  the  introduction  of  a trial  parameter  set  p0.  A new  set,  p, 
is  generated  via  equation  4.3  and  the  procedure  continues  until  the  error  function  is  below 
an  acceptable  tolerance. 

Performance  Of  Algorithms 

GAs  and  Newton  based  methods,  though  they  share  the  same  purpose,  are 
computationally  very  different.  Each  method  has  its  advantages  and  disadvantages  and 
as  a general  recommendation  for  the  success  of  the  optimization,  it  is  recommended  to 
use  a combination  of  the  two  algorithms. 

Perhaps  the  most  marked  advantage  of  the  GA  over  the  Newton  schemes  is  the 
global  minimum  convergence  of  the  GA  versus  the  Newton  algorithms.  Since  GAs  are 
stochastic  optimizations,  there  is  a random  “hopping”  around  the  surface  which  allows 
the  GA  to  converge  to  lower  minima  than  Newton  schemes  which,  given  an  initial  point, 
will  proceed  downhill  and  converge  locally.  Often,  Newton  schemes  get  trapped  in  local 
minima  and  have  no  provision  to  escape  these  areas  of  the  surface.  Although  GAs  may 
sample  these  local  minima  which  are  problematic  for  Newton  schemes,  due  to  the 
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random  nature  of  the  GA  search,  they  can  easily  escape  these  areas  of  the  surface  and 
find  better  solutions. 

Although  the  random  search  is  advantageous  as  demonstrated  in  the  previous 
discussion,  it  has  the  major  drawback  that  the  GA  requires  lengthy  amounts  of  computer 
time  to  converge.  A successful  GA  run  requires  a large  number  of  trial  parameter  sets 
from  which  to  select  fit  individuals.  Further,  the  GA  must  be  allowed  to  evolve  over 
many  generations  so  that  a suitable  set  of  parameters  can  be  found.  In  other  words,  in  the 
spirit  of  Darwin’s  theory  of  natural  selection,  evolution  must  be  allowed  to  proceed  over 
many  generations  in  order  to  find  the  most  fit  individual.  If  the  error  function  which  is 
being  minimized  can  be  rapidly  evaluated,  the  requirement  of  large  population  size  and 
multiple  generation  sampling  is  not  detrimental.  However,  in  the  parameterization  of  the 
transfer  Hamiltonian,  the  error  functions  are  computationally  costly,  relatively  speaking, 
requiring  the  evaluation  of  Hartree-Fock  energies  and  gradients  for  each  parameter  set 
which  results  in  a long  run  time  for  the  GA. 

One  final  deleterious  feature  of  GAs  is  that  after  many  generations,  if  the  global 
minimum  has  not  been  located,  the  search  will  begin  to  stagnate.  This  is  due  to  the 
successive  removal  of  comparatively  less  fit  individuals  from  the  population,  which,  after 
many  generations,  results  in  a population  of  parameter  sets  which  are  all  more  or  less 
equally  fit  and  the  evolution  of  the  parameters  stalls.  This,  in  a sense,  is  similar  to 
Newton  algorithms  getting  trapped  in  a local  minimum.  When  this  situation  occurs,  it  is 
advantageous  to  switch  algorithms  and  use  the  Newton  optimization  algorithm  to  further 
reduce  the  value  of  the  error  function.  This  combination  of  techniques  has  been 
successfully  applied  in  the  parameterizations  of  the  transfer  Hamiltonian. 
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Procedure 

It  is  difficult  to  give  a general  description  of  the  parameterization  procedure  for  the 
transfer  Hamiltonian  since  each  problem  behaves  differently  and  requires  problem 
specific  analysis.  In  general,  it  has  been  found  that  the  GA  should  be  initialized  with 
about  100  trial  parameter  sets  and  allowed  to  evolve  for  150  to  200  generations.  This 
gives  a balance  between  computational  time  and  sufficient  population  size  for  an 
effective  GA  run.  Using  the  resulting  converged  set,  the  parameters  should  be  refined 
using  the  Newton  type  algorithm  to  further  reduce  the  value  of  the  error  function.  For 
convergence  of  the  optimization,  it  is  always  best  to  focus  on  a smaller  set  of  reference 
data.  The  more  types  of  properties  targeted  in  the  parameterization,  the  worse  the 
transfer  Hamiltonian  will  perform  for  the  target  set  as  a whole.  A similar  conclusion  was 
reached  by  Bernstein  and  Kaxiras  in  their  attempt  to  parameterize  a tight  binding 
Hamiltonian  [31]. 

Transfer  Hamiltonian  Examples 

Since  the  goal  of  this  work  is  predictive  molecular  dynamics  simulation  which 
requires  accurate  interatomic  forces,  the  parameters  of  the  Hartree-Fock  variant  of  the 
transfer  Hamiltonian  operator  were  determined  such  that  the  interatomic  forces  at  the 
transfer  Hamiltonian  level  match  those  of  ab  initio  theory.  Specifically,  coupled  cluster 
(CCSD/DZP)  forces  were  determined  for  the  dissociation  reactions  of  several  small 
molecules  and  the  parameters  determined  such  that  the  reference  coupled  cluster  data  is 
reproduced  by  the  transfer  Hamiltonian. 

Ethane 

Coupled  cluster  forces  were  determined  for  the  dissociation  reaction  of  ethane. 


C2H6  «->  2CH3 


40 


where  the  carbon-carbon  bond  is  the  distinguished  coordinate  [58],  At  each  geometry, 
from  the  repulsive  region  out  to  dissociation,  the  C-C  bond  is  held  fixed  and  the 
remaining  3N-7  coordinates  allowed  to  relax.  The  resulting  force  curve  was  then  used  as 
reference  data  for  the  transfer  Hamiltonian. 

For  this  problem,  a minimal  basis  set  of  valence  (s,p)  orbitals  was  assigned  to 
each  carbon  atom  and  a single  s orbital  was  used  as  a basis  for  each  hydrogen.  In  this 
particular  optimization,  only  the  (a,b,c)  parameters  of  the  core  repulsion  function  (see 
chapter  3,  equation  3-21)  were  determined  in  the  parameterization  and  the  other 
parameters  were  held  fixed  at  their  given  MNDO/AM1  values  [59].  The  resulting 
parameters  for  the  core  repulsion  function  are  contained  in  the  appendix. 

Figure  4-1  is  a comparison  of  the  CCSD  force  curve,  versus  those  obtained  with 
the  transfer  Hamiltonian  and  standard  semi -empirical  theory.  In  the  figure,  it  is  seen  that 
the  transfer  Hamiltonian  parameterization  provides  an  accurate  representation  of  the  ab 
initio  force  curve. 

Disilane 

A similar  parameterization  was  performed  for  the  dissociation  of  disilane 

S^H^  *-*  2SiH3 

where  the  Si-Si  bond  is  the  distinguished  coordinate.  Similar  to  the  ethane  molecule, 
forces  were  generated  at  the  CCSD/DZP  level  of  theory  at  points  in  the  repulsive  and 
dissociation  regions  of  the  potential  energy  surface,  and  the  core  repulsion  parameters 
(AMI  Hamiltonian)  determined  such  that  the  transfer  Hamiltonian  forces  reproduce  those 
given  by  coupled  cluster  theory.  The  converged  parameters  are  contained  in  the  appendix 
and  the  resulting  force  curves  are  shown  in  Figure  4-2. 
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Comparison  Of  Forces  For  Ethane  Dissociation 
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Figure  4-1.  Comparison  of  forces  for  ethane  dissociation  using  coupled  cluster,  transfer 
Hamiltonian,  and  standard  semi-empirical  theory. 

As  figures  4-1  and  4-2  demonstrate,  the  parameterized  forces  for  the  transfer 
Hamiltonian  reproduce  those  given  by  the  ab  initio  CCSD  calculations.  It  should  also  be 
emphasized  that  the  forces  are  accurate  over  the  whole  extent  of  the  potential  energy 
surface,  and  not  just  near  equilibrium  as  most  of  the  existing  potentials  in  the  literature 
are.  A comparison  of  the  CPU  time  required  for  one  energy+gradient  calculation  for  the 
two  representative  systems  is  given  in  Table  4-1.  The  transfer  Hamiltonian  performs 
three  to  four  orders  of  magnitude  more  quickly  than  the  ab  initio  model,  yet  the  resulting 
forces  are  well  reproduced  as  seen  in  the  figures.  This  demonstrates  that  the  transfer 
Hamiltonian  provides  a convenient  route  for  the  inclusion  of  ab  initio  forces  in  molecular 


dynamics  simulations. 
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Comparison  Of  Forces  For  Disilane  Dissociation 
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Figure  4-2.  Comparison  of  forces  for  disilane  dissociation  using  coupled  cluster,  transfer 
Hamiltonian,  and  standard  semi -empirical  theory. 

Conclusion 

This  chapter  discussed  the  parameterization  strategy  of  the  transfer  Hamiltonian 
and  presented  two  test  systems  which  demonstrated  the  ability  of  the  model  to  reproduce 
ab  initio  forces,  over  a range  of  geometries,  at  a much  cheaper  computational  cost  than 
ab  initio  theory.  With  this  as  a framework,  the  remainder  of  the  dissertation  will  focus  on 
the  creation  of  transfer  Hamiltonians  for  silicon  oxide  systems  and  their  application  to  the 
fracture  of  a silica  model. 

Table  4-1.  Computer  time  (CPU  seconds)  necessary  to  perform  one  energy+gradient 
calculation  for  the  representative  systems. 


Molecule 

CCSD 

TH-CCSD 

c2h6 

45.1 

0.02 

Si2H6 

104.8 

0.02 

CHAPTER  5 

FRACTURE  SIMULATIONS 

Introduction 

This  chapter  will  discuss  the  parameterization  and  validation  of  a transfer 
Hamiltonian  applicable  to  the  polymorphs  of  silica.  It  will  be  seen  that  the  resulting 
transfer  Hamiltonian  performs  well  for  the  parameterization  cluster,  and  for  clusters 
outside  the  realm  of  the  parameterization.  This  parameterization  of  the  transfer 
Hamilitonian  will  then  be  used  in  a molecular  dynamics  simulation  of  fracture  in  a silica 
model  containing  108  quantum  mechanical  atoms. 

Fracture 

Fracture  begins  when  stress  is  applied  to  a material.  When  the  stress  exceeds  a 
critical  value,  atomic  bonds  begin  to  break  and  the  crack  propagates  in  the  system. 

The  failure  of  materials  due  to  fracture  is  a costly  and  undesirable  event  which  plagues 
many  industries[60].  Therefore,  an  understanding  of  the  underlying  physics  and 
chemistry  which  determine  the  macroscopic  failure  mechanisms  in  real  materials  is 
essential  so  that  premature  component  failure  can  be  predicted  and  alleviated. 

Fracture  may  be  classified  as  either  brittle  or  ductile.  In  general,  the  primary 
difference  between  the  two  fracture  modes  is  the  amount  of  plastic  deformation  (degree 
of  atomic  disorder)  that  the  material  undergoes  before  fracture  occurs.  Ductile  materials 
undergo  a large  degree  of  plasticity  whereas  brittle  materials  show  little  or  no  tendency  to 
do  so.  In  ductile  materials,  once  the  crack  is  initiated,  it  moves  slowly  and  will  usually 
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not  extend  unless  an  increased  stress  is  applied.  On  the  contrary,  in  brittle  fracture, 
cracks  spread  very  rapidly  and  are  accompanied  by  a low  energy  release  [61]. 

Typically,  it  is  regarded  that  ductile  fracture  is  more  desirable  than  brittle  fracture 
in  terms  of  structural  or  component  failure  [62].  Since  ductile  fracture  occurs  over  a 
period  of  time  and  is  generally  a slower  process  than  brittle  fracture,  possible  sites  of 
material  fatigue  can  be  identified  and  repaired  before  catastrophic  failure. 

Transfer  Hamiltonian  - Silica 

Silica  is  a prototypical  brittle  material.  Given  the  importance  of  silica  from  a 
geological  and  technical  standpoint,  the  identification  of  the  mechanical  properties  of  this 
material  has  resulted  in  a significant  number  of  molecular  dynamics  studies  [63-66].  In 
many  of  these  studies,  the  BKS  or  TTAM  potentials  discussed  in  the  second  chapter  of 
this  thesis  were  used  to  describe  the  interatomic  forces  which  are  the  determining  factor 
for  a predictive  simulation.  However,  as  demonstated  previously,  these  potentials  fail  to 
reproduce  the  forces  obtained  from  high  level  ab  initio  calculations,  therefore  results 
based  on  these  models  are  questionable.  On  the  contrary,  as  demonstrated  in  Chapter  4 
for  ethane  and  disilane,  the  transfer  Hamiltonian  reproduces  coupled  cluster  theory  forces 
at  a fraction  of  the  computational  cost  of  the  ab  initio  methods  which  makes  it  a useful 
tool  for  molecular  dynamics  and  enables  a predictive  simulation  of  silica  to  be  performed, 
given  a suitable  parameterization.  What  remains  then  is  a generation  of  ab  initio 
reference  data  for  some  small  model  of  silica,  the  parameterization  of  the  transfer 
Hamiltonian  based  on  the  reference  data,  and  finally,  the  actual  simulation  of  a much 
larger  sized  system  using  the  transfer  Hamiltonian  as  a potential  for  the  dynamics. 
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Reference  Cluster 

In  the  parameterization  of  any  potential,  a suitable  set  of  reference  data  has  to  be 
generated  and  the  potential  parameters  fit  such  that  the  reference  data  is  reproduced  by 
the  potential.  In  practice,  for  molecular  dynamics  simulations  of  systems  which  contain  a 
large  number  of  atoms,  a small  model  of  the  target  simulation  system  is  selected  and  the 
relevant  properties  of  the  smaller  molecule  used  as  a reference  data  set.  After  the 
parameterization,  the  resulting  potential  is  then  assumed  to  be  valid  for  systems  outside 
the  reference  data  set,  such  as  the  target  simulation  system.  Given  that  assumption, 
simulations  of  very  large  systems  can  performed  using  that  potential  to  compute 
interatomic  forces. 

In  this  procedure,  it  is  clear  that  the  choice  of  parameterization  cluster  is 
extremely  important.  In  general,  the  cluster  should  provide  a good  representation  of  the 
local  geometry  of  the  target  simulation  system.  Silica  is  known  to  assume  a wide  variety 
of  crystal  structures  such  as  quartz,  cristobalite,  coesite,  and  stishovite.  Although  the 
density  of  each  polymorph  varies,  they  all  have  a region  of  stability,  given  the  proper 
ambient  conditions.  A common  feature  of  all  the  stable  structures  is  the  tetrahedral  SiC>4 
building  unit  where  a central  silicon  is  surrounded  by  four  oxygen  atoms.  The  tetrahedra 
share  one  or  more  oxygens  with  neighboring  tetrahedra  and  their  relative  orientation 
defines  the  type  of  structure.  This  tetrahedral  type  of  bonding  is  also  present  in 
amorphous  silica.  Given  the  importance  of  the  tetrahedral  building  block  in  the 
polymorphs  of  silica,  pyrosilicic  acid,  Si207H6  (see  Figure  5-1),  has  been  chosen  as  a 
small  molecular  model  of  bulk  silica.  This  molecule  has  been  used  in  cluster  studies  of 
silica  by  other  researchers  as  well  [67-68], 
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Figure  5-1.  Structure  of  pyrosilicic  acid. 

Reference  Data  For  Pyrosilicic  Acid 

Given  the  pyrosilicic  acid  molecule,  ab  initio  coupled  cluster  (CCSD/DZP) 
calculations  were  performed  to  generate  the  reference  data  to  which  the  transfer 
Hamiltonian  parameters  were  fit.  For  the  coupled  cluster  computations,  the  p orbital 
polarization  functions  of  the  DZP  basis  were  removed  from  the  hydrogens  and  the  Si-O- 
Si  bond  angle  was  fixed  at  180°.  The  latter  restriction  was  enforced  to  maintain 
symmetry  and  reduce  computational  effort  for  the  coupled  cluster  calculations.  Initially, 
the  entire  system  was  allowed  to  relax  under  the  imposed  symmetry  constraints  until  an 
equilibrium  configuration  was  found.  Starting  from  the  relaxed  geometry,  one  of  the 
equivalent  bridged  Si-0  bonds  was  chosen  as  the  distinguished  coordinate  and  the  forces, 
both  on  the  compressive  and  dissociative  sides  of  the  equilibrium  bond  length,  were 
determined.  For  the  determination  of  the  forces  along  the  distinguished  coordinate  path, 
the  chosen  Si-0  bond  is  held  fixed  at  varying  lengths  and  the  other  internal  degrees  of 
freedom,  except  those  of  the  hydrogen  atoms,  allowed  to  relax.  The  constraint  on  the 
hydrogens  was  enforced  to  prevent  any  spurious  hydrogen  bonding  interactions  which 
would  not  be  present  in  real  silica.  The  resulting  forces  from  these  constrained 
optimizations  are  used  as  reference  data  for  the  transfer  Hamiltonian.  Given  an  accurate 
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fit  for  the  ab  initio  forces,  the  transfer  Hamiltonian  can  then  be  applied  to  large  molecular 
dynamics  simulations  of  silica. 

Parameterization 

For  this  study,  the  Hartree-Fock  variant  of  the  transfer  Hamiltonian  was  used.  A 
minimal  valence  basis  set  (s,p)  was  assigned  to  silicon  and  oxygen  in  pyrosilicic  acid  and 
no  polarization  functions  were  assigned  to  the  hydrogen  atoms.  The  one  and  two  electron 
integrals  were  computed  using  the  formulas  presented  in  Chapter  3 with  the  parameters 
held  fixed  at  their  existing  values  (AMI  parameters)  [59,69].  Only  the  (a,b,c)  parameters 
of  the  core  repulsion  function  (see  equation  3-21)  were  varied  in  the  optimization.  The 
optimized  parameters  of  the  silicon  and  oxygen  atoms  in  pyrosilicic  acid  were  determined 
using  the  genetic  and  Newton  based  algorithms  such  that  the  coupled  cluster  forces  along 
the  same  distinguished  coordinate  path  used  in  the  ab  initio  work  was  accurately 
reproduced  by  the  transfer  Hamiltonian.  The  resulting  parameter  set  is  contained  in  the 
appendix. 

Validation  Of  Potential 

Given  the  converged  parameter  set,  the  fidelity  of  the  potential  must  be  examined 
to  ensure  that  it  is  indeed  reliable.  The  potential  must  be  able  to  properly  reproduce  the 
given  reference  data  and  should  also  be  able  to  describe  systems  which  were  not 
explicitly  considered  in  the  parameterization.  Several  tests  of  the  current 
parameterization  of  the  transfer  Hamiltonian  for  silica  were  conducted  and  are  presented 
in  the  following. 

Figure  5-2  is  a plot  of  the  forces  obtained  for  pyrosilicic  acid  dissociating  into 
neutral  fragments  via 


Si207H6  «-►  OSi(OHh  + Si(OH)3 


(5.1) 
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using  coupled  cluster,  transfer  Hamiltonian  and  BKS  potentials.  In  the  figure  it  is  seen 
that  the  transfer  Hamiltonian  reproduces  the  coupled  cluster  forces  in  the  repulsive  and 
dissociative  regions  of  the  potential  energy  surface,  unlike  the  classical  potential.  This  is 
an  important  feature  since  these  potentials  are  used  to  study  fracture  where  interatomic 
bonds  are  stretched  from  their  equilibrium  values  and  need  proper  description  over  the 
entire  range  of  bond  lengths. 


Figure  5-2.  Comparison  of  forces  for  dissociation  of  pyrosilicic  acid  into  neutral 
fragments  using  ab  initio  and  parameterized  potentials. 

The  forces  plotted  in  figure  5-2  pertain  to  the  lower  energy  dissociation  pathway 
of  pyrosilicic  acid  into  neutral  fragments  and  the  reference  data  to  which  the  transfer 
Hamiltonian  parameters  were  fit  was  determined  for  this  particular  dissociation  path. 
However,  the  system  may  also  separate  into  charged  fragments  via 

Si207H6*->  OSi(OH)3  + +Si(OH)3  (5.2) 

where  the  different  pathways  can  be  studied  computationally  by  using  a restricted 
(charged)  or  unrestricted  (neutral)  Hartree-Fock  computation  to  generate  the  ground  state 
reference  determinant.  Dissociation  into  charged  species  lies  ==30  kcal/mol  higher  in 
energy  than  the  neutral  dissociation  and  as  a test  on  the  veracity  of  the  parameter  set 
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obtained  for  this  model  cluster,  the  transfer  Hamiltonian  forces  for  the  charged 
dissociation  were  compared  to  those  computed  with  coupled  cluster  theory.  A plot  of  the 
forces  for  the  dissociation  into  charged  fragments  is  shown  in  Figure  5-3.  There  it  is  seen 
that  the  coupled  cluster  forces  for  dissociation  into  charged  fragments  are  recovered  by 
the  transfer  Hamiltonian,  although  the  parameterization  only  included  data  from  the 
lower  energy  pathway.  This  type  of  flexibility  is  not  possible  with  a classical  potential 
such  as  BKS  or  TTAM  which  does  not  recognize  the  difference  between  ions  and  neutral 
species. 


Figure  5-3.  Comparison  of  forces  for  pyrosilicic  acid  dissociation  into  charged  fragments. 

The  potential  energy  surface  of  the  two  dissociation  pathways  is  shown  in  Figure  5- 
4.  The  energy  differences  for  the  transfer  Hamiltonian  were  obtained  by  integration  of 
the  forces  computed  with  the  converged  parameter  set.  As  the  figure  demonstrates,  the 
potential  energy  surface  along  each  pathway  is  reproduced  by  the  transfer  Hamiltonian 
although  no  explicit  consideration  of  energies  was  included  in  the  parameterization. 
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Figure  5-4.  Comparison  of  potential  energy  surfaces  for  pyrosilicic  acid  dissociation 
into  neutral  and  charged  fragments  using  coupled  cluster,  transfer 
Hamiltonian  and  standard  semi-empirical  theory. 
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In  order  for  the  transfer  Hamiltonian  to  be  applicable  to  large  polymorphs  of  silica, 
it  is  necessary  that  the  parameters  be  convergent  with  respect  to  cluster  choice.  In  other 
words,  the  parameters  which  were  determined  for  pyrosilicic  acid  must  not  change 
(significantly)  if  the  next  larger  sized  Si-0  cluster,  trisilicic  acid,  is  used  for  the 
parameterization.  If  the  parameters  are  indeed  convergent,  then  it  can  be  argued  that  the 
potential  is  applicable  to  larger  Si-0  systems  of  all  sizes.  This  is  achievable  in  the 
context  of  the  transfer  Hamiltonian  because  the  operator  is  very  short  ranged  (limited  to 
two  center  interactions  at  most)  and  the  interatomic  interactions  fall  off  very  rapidly  with 
distance.  Therefore,  the  inclusion  of  additional  Si-0  units  which  are  long  distances  from 
the  site  of  interest  have  a nominal  effect  on  the  local  forces.  To  exemplify  this,  consider 
Figure  5-5  which  is  a comparison  of  the  forces  for  trisilicic  acid,  which  contains  one 
more  Si04  unit  than  does  the  pyrosilicic  acid  parameterization  cluster.  Due  to  the  larger 
cluster  size,  this  system  was  treated  with  DFT,  however,  the  resulting  potential  energy 
surface  and  forces  should  closely  follow  those  that  would  be  obtained  with  coupled 
cluster  theory.  As  shown,  the  transfer  Hamiltonian  force  curve  generated  with 
parameters  determined  for  a smaller  cluster,  is  a good  representation  of  the  larger  cluster 
forces.  This  would  also  apply  for  even  larger  clusters  of  SiC>4  units  due  to  the  local 
nature  of  the  transfer  Hamiltonian  and  is  evidence  that  this  potential  is  applicable  to  bulk 
silica  samples. 

The  CPU  time  required  for  each  energy+gradient  evaluation  using  pyrosilicic  acid 
is  shown  in  Table  5-1.  The  computational  time  for  the  ab  initio  coupled  cluster 
calculation  is  4 orders  of  magnitude  more  costly  than  the  transfer  Hamiltonian 
computation  although  the  forces  and  energy  differences  are  equivalent  in  both  models. 


52 


Figure  5-5.  Comparison  of  forces  for  trisilicic  acid. 

Although  the  transfer  Hamiltonian  performs  quickly  enough  for  use  in  simulations, 
as  shown  in  the  table,  the  classical  BKS  potential  is  still  much  faster  than  the  transfer 
Hamiltonian.  The  simplicity  of  the  classical  potential  form  renders  it  extremely  fast  and 
the  transfer  Hamiltonian,  although  a step  in  the  right  direction,  will  never  be  able  to 
compete  with  the  must  simpler  classical  potential  forms  in  regard  to  computational  speed. 
Table  5-1.  Computer  time  (CPU  seconds)  necessary  to  perform  one  energy+gradient 


calculation  for  the  pyrosilicic  acic 

Method 

system. 

CPU  Time  (s) 

CCSD 

8656 

DFT 

375 

Transfer  Hamiltonian 

0.17 

BKS 

0.001 

As  stated  previously,  for  potentials  of  both  quantum  and  classical  forms,  it  is 
assumed  that  the  converged  parameter  set  is  valid  for  systems  outside  the  realm  of  the 
reference  data  set.  This  assumption  has  to  be  made  because  it  is  obviously  impossible  to 
include  every  possible  situation  in  one  parameterization.  Although  this  “leap  of  faith”  is 
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unavoidable,  it  is  still  important  to  assess  the  applicability  of  the  potentials  to  systems 
which  differ  from  that  to  which  they  were  parameterized.  A faithful  potential  should  be 
able  to  handle  a variety  of  bonding  situations,  especially  potentials  which  are  used  to 
study  silica.  The  surface  of  amorphous  silica  is  known  to  contain  many  defect  sites 
which  includes  under  coordinated  silicon  atoms  and  small,  highly  strained  rings.  A 
predictive  simulation  of  silica  must  use  a potential  which  is  capable  of  properly 
characterizing  all  these  bonding  patterns. 

Figure  5-6  contains  histograms  detailing  the  error  between  bond  lengths  and  angles 
in  several  small  silicon  oxide  systems  different  from  pyrosilicic  acid.  In  comparison  to 
CCSD(T)  geometries,  the  transfer  Hamiltonian  performs  very  well  for  these  systems, 
although  none,  some  of  which  contain  double  bonds,  were  included  in  the 
parameterization  set.  Also  shown  in  the  figure  is  the  geometric  error  obtained  with  the 
TTAM  and  LSU  3-body  potential.  The  error  in  these  internal  coordinates  using  these 
potentials  is  greater  than  that  with  the  transfer  Hamiltonian,  and  the  LSU  potential  fails  to 
produce  the  linear  bond  arrangement  in  the  triatomic  molecule  SiC>2.  This  is  due  to  the 
preferential  treatment  of  tetrahedral  bond  angles  built  into  the  LSU  potential  (see  Chapter 
2). 

As  a final  test  on  the  robustness  of  the  parameterization,  the  transfer  Hamiltonian 
was  applied  to  a class  of  systems  known  as  polyhedral  oligomeric  silsesquioxanes 
POSS).  These  cage-like  silicon-oxide  molecules  are  closely  related  to  silica  and  silicones 
and  have  shown  great  promise  as  additives  which  greatly  improve  the  thermal  and 
physical  properties  of  many  plastics  [70].  Xiang  et  al.  [71]  performed  a DFT  study 
(of  several  of  these  systems,  and  presented  the  structure  and  energy  ordering  for  POSS 
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Comparison  of  Si-0  Bond  Lengths  Relative  to 
CCSD(T) 


Comparison  of  O-Si-O  Angle 


Figure  5-6.  Comparison  of  bond  lengths  and  angles  for  several  small  silicon  oxide 
systems  using  various  potentials. 


55 


isomers  with  equivalent  stoichiometries.  The  eight  and  twelve  silicon  atom  cages  shown 
in  Figure  5-7  were  treated  with  the  transfer  Hamiltonian  to  determine  if  the  energy 
ordering  among  the  isomers  given  by  DFT  was  reproduced  by  the  transfer  Hamiltonian. 


Figure  5-7.  Structure  of  the  eight  and  twelve  member  POSS  cages.  The  isomers  are 

labeled  by  their  point  group  symmetries  and  number  of  silicon  atoms  in  each 
cage. 

Xiang  et  al.  report  that  among  the  two  eight  member  isomers,  the  conformer  of  Oh 
symmetry  is  the  most  stable.  For  the  twelve  member  cage,  the  D2d  isomer  is  the  most 
favorable  geometry.  The  transfer  Hamiltonian  reproduces  these  structures  and  energetic 
orderings. 

Given  the  performance  of  the  transfer  Hamiltonian  in  the  tests  given  above,  it  can 
be  concluded  that  the  transfer  Hamiltonian  is  a useful  tool  for  the  study  of  silica.  The 
forces  for  pyrosilicic  acid  match  those  of  coupled  cluster  theory  in  the  repulsive  and  bond 
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break  regions  which  is  necessary  for  fracture  simulations.  Further,  the  description  of 
other  silicon  oxide  systems  of  different  coordinations  is  well  described  which  is 
necessary  to  properly  model  defect  sites  in  silica. 

Application 

The  transfer  Hamiltonian  was  used  to  study  the  failure  response  to  uniaxial  tension 
of  a 108  atom  silicon-oxide  nanorod  (Figure  5-8)  via  molecular  dynamics  simulation. 
The  nanorod  is  a model  system  of  silica  which  is  large  enough  to  yield  genuine  materials 
properties  such  as  stress-strain  curves  yet  can  be  treated  fully  quantum  mechanically  at 
the  transfer  Hamiltonian  level. 

The  nanorod,  as  shown  in  the  figure,  consists  of  planar  six  membered  rings, 

Each  planar  ring  shares,  both  above  and  below,  a ring  of  oxygen  atoms  with  adjacent  six 
member  rings.  To  alleviate  any  dangling  bonds,  the  nanorod  is  terminated  with  capped 
rings  where  each  silicon  atom  of  the  terminating  ring  is  connected  by  a bridging  oxygen 
or  two  interstitial  oxygens  as  shown.  The  nanorod  has  the  proper  stoichiometric  ratio  of 
silicon  to  oxygen  observed  in  real  silica  (1:2)  and  is  considered  to  be  a viable  model  for 
the  study  of  silica.  In  this  particular  work,  the  nanorod  was  used  to  elucidate  the 
differences,  if  any,  of  the  fracture  mode  of  the  nanorod  using  classical  versus  quantum 
potentials. 

Fracture  Of  Nanorod 

Before  presentation  of  the  simulations  and  as  an  additional  demonstration  of  the 
faithful  applicability  of  the  transfer  Hamiltonian  to  large  systems,  the  energy  of  the 
nanorod  was  computed  at  equilbrium  and  strained  geometries  using  the  transfer 
Hamiltonian  and  density  functional  theory.  The  nanorod  was  equilibrated  with  each 
potential  and  successively  expanded  by  a few  percent  to  generate  the  energy  profiles 
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Figure  5-8.  108  atom  silica  nanorod  viewed  from  side  and  top.  The  larger  circles  are 
silicon  atoms  and  the  others  are  oxygen. 

shown  in  Figure  5-9.  The  potential  curves  show  fairly  good  agreement  considering  the 
initial  parameterization  focused  on  accurate  forces  for  a much  smaller  cluster.  Although 
the  agreement  with  DFT  is  good,  the  transfer  Hamiltonian  operates  much  more  quickly. 
Each  single  point  computation  required  43  CPU  seconds  for  the  transfer  Hamiltonian, 
while  the  DFT  computations  used  85,000  seconds  of  CPU  time. 

To  study  fracture,  the  nanorod  was  subjected  to  uniaxial  tension  using  molecular 
dynamics  simulation  where  the  stress  was  applied  by  fixing  the  longitudinal  movement  of 
the  fifteen  atoms  contained  in  the  caps  on  both  ends  of  the  nanorod.  With  the  long  axis 
of  the  nanorod  considered  to  be  the  z axis,  each  end  cap  atom  was  given  a fixed  velocity 
in  the  loading  direction  (along  z axis)  and  the  remaining  atoms  allowed  to  relax,  with 
their  positions  computed  by  the  deterministic  molecular  dynamics  algorithms.  The 
simulations  were  conducted  at  a constant  temperature  of  IK  using  the  predictor-corrector 
algorithm,  where  the  temperature  was  maintained  by  velocity  rescaling  [3,4]. 
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Figure  5-9.  Plot  of  potential  energy  curves  for  the  nanorod  using  transfer  Hamiltonian 
and  density  functional  potentials. 

The  stress-strain  curve  obtained  for  fracture  of  the  nanorod  using  the  transfer 
Hamiltonian  is  shown  in  Figure  5-10.  Curves  obtained  using  BKS  and  TTAM  classical 
potentials  are  included  for  comparison.  The  stresses  plotted  in  the  figure  correspond  to 
the  sum  of  all  forces  parallel  to  the  loading  direction  in  the  constrained  atoms  of  the  end 
caps  divided  by  the  projected  cross  sectional  area  of  the  end  caps.  The  cross  sectional 
area  is  updated  at  each  time  step  and  accounts  for  any  effects  due  the  xy  motion  of  the 
atoms  in  the  longitudinally  constrained  end  caps.  This  differs  from  the  usual  engineering 
stress-strain  curve  where  the  initial  cross  sectional  area  is  used  at  each  time  step.  The 
strain  is  computed  as  (L-L0)/L0,  where  L and  L0  are  the  current  and  initial  lengths  of  the 
nanorod  respectively. 

Examination  of  the  stress-strain  curves  indicates  that  the  transfer  Hamiltonian  and 
classical  potentials  predict  approximately  equivalent  values  for  the  critical  strain, 
however  the  tensile  strengths  for  the  transfer  Hamiltonian  and  BKS  potentials  are  higher 
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than  those  predicted  with  the  TTAM  potential.  The  curvature  of  the  classical  potential 
curves  indicates  that  the  nanorod  undergoes  a more  plastic  deformation,  characteristic  of 
ductile  fracture,  than  that  of  the  transfer  Hamiltonian  which  maintains  an  essentially 
linear  stress-strain  relationship  (brittle  fracture)  up  to  failure. 

The  qualitatively  distinct  fracture  modes  indicated  by  the  classical  (ductile)  and 
quantum  potential  (brittle)  exemplifies  the  importance  of  having  accurate  forces  in  the 
dynamics  simulation.  It  cannot  be  concluded  which  failure  mode  is  actually  correct,  but 
given  the  performance  of  the  transfer  Hamiltonian  for  the  tests  discussed  previously,  the 
results  it  yields  for  the  nanorod  should  be  considered  the  most  reliable.  The  failure  of  the 
classical  potentials  to  properly  characterize  the  mechanical  behavior  of  the  nanorod 
could,  were  the  nanorod  of  technological  importance,  result  in  wasted  time  and  effort  for 
experimentalists  due  to  misleading  recommendations  by  simulators.  Of  course  this  is  a 
hypothetical  situation  at  present,  however,  the  implications  are  clear. 

Hybrid  Dynamics  Of  Nanorod 

The  CPU  time  required  for  the  nanorod  simulations  is  presented  in  Table  5-2.  The 
quantum  simulations  require  much  more  time  to  reach  completion  than  do  the  classical 
potentials.  This  difference  will  become  even  more  pronounced  as  the  size  of  the 
simulation  system  increases  and  the  quantum  calculation  becomes  swamped  by  the  huge 
dimension  of  the  matrices  involved  in  the  calculation.  This  situation  can  be  rectified  by 
the  use  of  a hybrid  molecular  dynamics  scheme  whereby  the  forces  of  interaction  are 
computed  using  the  much  more  rapid  classical  potential  and  the  slower  transfer 


Hamiltonian  in  the  same  simulation. 
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Figure  5-10.  Stress-Strain  curves  for  uniaxial  tension  of  108  atom  nanorod. 
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Table  5-2.  Computer  time  (CPU  seconds)  necessary  to  perform  5000  timestep  MD 


simulation. 

Potential 

CPU  Time 

Transfer  Hamiltonian 

2.5  Days 

TTAM/BKS 

1.5  Minutes 

In  the  hybrid  simulation  scheme,  the  transfer  Hamiltonian  is  used  to  compute  the 
force  vector,  Fth,  at  the  current  configuration.  Those  forces  are  stored  and  the  classical 
potential  force  vector,  Fcl,  at  the  same  configuration  is  determined.  The  difference 
between  the  two,  Ferr=Fth-Fcl,  is  a vector  of  errors  between  the  quantum  forces  and  the 
classical  forces.  Given  these  errors,  the  classical  potential  is  used  to  generate  forces  for 
the  next  few  configurations  and  the  previously  determined  quantum  corrections  are  added 
to  each  component  of  the  classical  force  vector  at  each  time  step.  This  scheme  has  the 
advantage  that  fewer  of  the  more  costly  transfer  Hamiltonian  calculations  need  to  be 
performed  since  the  majority  of  the  time  steps  are  evaluated  with  the  much  faster 
classical  potential. 

The  accuracy  of  this  force  correction  scheme  was  assessed  by  simulating  the 
fracture  of  the  nanorod,  as  before,  but  using  the  classical  TTAM  potential  with  periodic 
force  corrections  from  the  transfer  Hamiltonian.  Three  simulations  were  performed 
where  the  transfer  Hamiltonian  force  corrections  were  evaluated  at  intervals  of  10,  25, 
and  50  time  steps.  The  resulting  stress-strain  curves  are  shown  in  Figures  5-11  through 
5-13,  as  a function  of  the  frequency  of  quantum  corrections  performed  in  the  simulation. 
As  the  figures  show,  the  force-corrected  TTAM  results  are  closer  to  the  pure  transfer 
Hamiltonian  curve,  although  the  scheme  does  introduce  a large  amount  of  numerical 
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“noise.”  As  the  frequency  of  quantum  corrections  increases,  the  force-corrected  TTAM 
simulation  approaches  that  of  the  transfer  Hamiltonian  more  closely,  the  limit  being  a 
quantum  correction  at  each  timestep,  i.e.  a fully  quantum  simulation. 

Given  the  decreased  number  of  time  steps  for  which  the  transfer  Hamiltonian 
potential  must  be  evaluated,  the  hybrid  simulation  gives  a dramatic  speedup  in 
computational  time,  as  shown  in  Table  5-3  which  lists  the  CPU  time  required  for  each 
simulation. 

Conclusion 

This  chapter  discussed  the  parameterization  and  validation  of  a transfer  Hamiltonian 
applicable  to  silica.  It  was  shown  that  the  potential  is  accurate  for  the  parameterization 
cluster  and  other  silicon-oxide  systems  outside  the  parameterization  set,  including  POSS 
systems.  The  fracture  mode  of  a 108  atom  silicon  oxide  nanorod  was  simulated  and  the 
transfer  Hamiltonian  yielded  a qualitatively  different  material  behavior  than  did  the 
classical  BKS  and  TTAM  potentials.  Finally,  a computationally  efficient  hybrid 
simulation  scheme  was  presented  wherein  the  classical  force  fields  are  updated  with 
quantum  information  which  yields  a decreased  simulation  time  and  congruent  stress- 


strain  curves. 

Table  5-3.  CPU  time  for  hybrid  simulations. 

Force  Correction  Interval  (Timesteps) 

CPU  Time 

0 (Fully  Quantum) 

2.5  Days 

10 

1 .6  Hours 

25 

1 Hour 

50 

1/2  Hour 

250 
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Figure  5-11.  Stress  strain  curves  using  a quantum  force  correction  every  50  timesteps. 


250 


CO 

o 


co 

o 


Lfi 

CVJ 

O 


CVJ 

O 


LD 

O 


in 

o 


o 


o 

o 

CvJ 


o 

Ln 


o 

o 


o 

in 


o 


o 

in 

i 


c 

'E 

-4-» 

CO 


(Ed0)  SS&J4S 


Figure  5-12.  Stress  strain  curves  using  a quantum  force  correction  every  25  timesteps. 
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Figure  5-13.  Stress  strain  curve  using  quantum  correction  every  10  timesteps. 


CHAPTER  6 

HYDROLYTIC  WEAKENING 

Introduction 

This  chapter  will  discuss  the  changes  in  the  mechanical  properties  of  silica  in  the 
presence  of  water  by  a phenomena  known  as  hydrolytic  weakening.  A recently  presented 
mechanism  for  this  process  is  outlined  and  encoded  into  the  parameters  of  a transfer 
Hamiltonian.  The  transfer  Hamiltonian  will  then  be  used  in  molecular  dynamics 
simulations  of  water  assisted  fracture  of  the  nanorod  and  silicon-oxide  rings. 

Hydrolytic  Weakening 

It  is  well  known  that  water  greatly  enhances  the  rate  of  crack  growth  in  silica  in  a 
process  known  as  hydrolytic  weakening  [12].  Although  the  mechanism  of  this  process  is 
still  under  debate,  it  is  generally  accepted  that  this  is  a stress  induced  phenomenon 
whereby  water  assists  the  rupture  of  Si-0  bonds  in  the  silica  structure  and  enhances  the 
fracture  process.  One  model  of  this  process,  due  to  Michalske  and  Freiman  [73], 
identifies  the  water  monomer  as  the  primary  mediator  of  water  assisted  rupture  of  silicon- 
oxygen  bonds  in  silica.  In  this  mechanism,  it  is  suggested  that  a single  water  molecule 
orients  a lone  pair  of  electrons  on  its  oxygen  towards  a silicon  atom  of  the  surface.  A 
subsequent  proton  transfer  from  water  to  an  oxygen  of  the  silica  surface  results  in  rupture 
of  the  surface  Si-0  bond  and  enhanced  crack  tip  propagation. 

The  Michalske-Freiman  model  has  been  examined  by  Lindsay  et  al.  [68]  who 
applied  Hartree-Fock  level  calculations  to  pyrosilicic  acid  in  the  presence  of  several 
molecules,  including  water  and  ammonia.  Their  work  indicates  that  the  approach  of  the 
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monomer,  in  the  manner  described  by  Michalske  and  Freiman,  is  favorable,  however  the 
conclusions  given  in  the  paper  are  limited  by  the  computational  methodology  applied  to 
the  test  system. 

To  further  elucidate  the  water  assisted  bond  rupture  process,  Del  Bene  et  al.  [74] 
applied  large  basis  set,  correlated  methods  to  H3SiOSiH3  in  the  presence  of  the  water 
monomer,  dimer,  and  trimer.  In  this  work,  it  was  found  that  the  binding  energies  of  the 
water  clusters  with  the  H3SiOSiH3  increased  in  the  order  monomer<trimer<dimer  which 
indicates  that  the  water  dimer  interacts  more  favorably  with  the  Si-0  bond  than  does  the 
monomer  used  in  the  Michalske-Freiman  model.  Del  Bene  et  al.  also  identified  a 
mechanism  of  Si-0  bond  rupture  using  the  water  dimer.  In  this  scheme,  two  of  the 
protons  in  H3SiOSiH3  were  replaced  by  fluorine  atoms  (see  Figure  6-1)  to  model  the 
effect  of  the  bulk  in  large  silica  samples.  As  the  Si-O  bond  was  stretched,  a dual  proton 
transfer  occurred  which  mediated  the  rupture  of  the  bond  as  outlined  in  the  figure. 

Given  this  mechanism  for  water  assisted  bond  cleavage  in  silica,  what  remains  is  the 
encoding  of  this  information  into  a transfer  Hamiltonian  parameter  set  and  the 
demonstration  of  the  process  by  molecular  dynamics  simulation. 


Figure  6-1.  Schematic  representation  of  dual  proton  transfer  mechanism. 
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Transfer  Hamiltonian 

Parameterization 

To  successfully  simulate  the  interaction  of  water  with  silica  given  the  mechanism 
described  above,  the  forces  which  describe  water-water,  water-silica,  and  silica-silica 
interactions  have  to  be  described  quantitatively.  The  transfer  Hamiltonian  described  in 
Chapter  5 only  characterizes  the  silica-silica  interactions  and  does  not  accurately 
reproduce  the  interactions  of  silica  with  water,  therefore  a reparameterization  is 
necessary. 

The  parameterization  of  this  mechanism  is  very  difficult  given  the  requirements 
imposed  upon  the  transfer  Hamiltonian  parameter  set.  In  this  case,  as  for  the  previous 
systems,  accurate  molecular  dynamics  is  the  goal  therefore  the  parameterization  will 
focus  on  accurate  reproduction  of  ab  initio  forces  for  dissociation  of  water,  and  silica  in 
the  presence  of  water.  The  water  molecules  are  composed  of  oxygen/hydrogen  atoms 
and  silica  contains  silicon/oxygen.  Therefore,  there  are  three  atomic  species  which  are 
available  for  parameterization.  This  parameterization  can  proceed  in  two  ways.  First,  the 
entire  reference  data  set  (forces)  for  the  water-water  and  water-silica  interactions  can  be 
targeted  in  one  huge  parameterization  run  which  includes  water  and  the  water-silica 
cluster  forces  at  several  geometries.  An  alternative  is  to  parameterize  the  water  cluster 
interaction  separately,  and  then,  with  the  fixed  oxygen/hydrogen  parameters,  determine 
the  silicon  parameters  for  the  silica  model  cluster.  The  latter  option  is  favorable  because 
although  two  successive  optimizations  must  be  performed,  the  target  data  set  is  smaller 
for  each  run  and  yields  improved  convergence. 

Proceeding  in  this  manner,  converged  oxygen/hydrogen  parameters  were 
determined  which  reproduce  coupled  cluster  (CCSD)  geometries,  dipole  moments,  and 
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forces  for  the  water  monomer  and  dimer  using  the  Hartree-Fock  variant  of  the  transfer 
Hamiltonian.  The  MBPT(2)  forces  for  the  silica  model  (Figure  6-1)  were  then  fit  by 
variation  of  the  silicon  parameters  using  fixed  values  for  the  oxygen/hydrogen 
parameters.  The  transfer  Hamiltonian  description  of  oxygen  and  silicon  includes  addition 
of  spherical  harmonic  d functions  and  a set  of  p functions  for  each  hydrogen  in  addition 
to  the  usual  minimal  basis,  valence  set  of  basis  functions.  The  AMI  core  repulsion 
function  was  also  included  to  increase  the  flexibility  of  the  fit  in  the  core  terms.  The 
converged  parameter  set  is  contained  in  the  appendix. 

Validation 

Selected  internal  coordinates  and  the  dipole  moments  for  the  water  monomer/dimer 
are  contained  in  Table  6-1  where  it  is  seen  that  the  CCSD  geometry  of  the  monomer  and 
dimer  is  well  reproduced  by  the  transfer  Hamiltonian.  The  computed  dipole  moment  for 
the  monomer  is  accurate  but  that  for  the  dimer  is  in  error  by  0.2D. 

Table  6-1.  Selected  internal  coordinates  for  the  water  monomer  and  dimer.  The  CCSD 


values  are  given  in  parenthesis.  All  values  in  Angstroms,  degrees,  and  Debye. 


Coordinate 

Monomer 

Dimer 

R(OH) 

0.96(0.95) 

0.96(0.95) 

0(HOH) 

105(104) 

105(102) 

0(HOH)b 

100(101) 

R(OO) 

2.90(2.90) 

Dipole 

1.85(1.86) 

2.45(2.68) 

Figure  6-2  is  a plot  of  the  forces  for  dissociation  of  the  water  monomer  using  the  transfer 
Hamiltonian  and  coupled  cluster  theory  for  the  reaction 


H2O^OH+H 


(6.1) 
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The  transfer  Hamiltonian  accurately  reproduces  the  coupled  cluster  forces  in  both  the 
repulsive  and  dissociative  regions  of  the  potential  energy  surface  which  is  required  for 
the  proton  transfer  to  be  modeled  accurately. 


Figure  6-2.  Comparison  of  forces  for  dissociation  of  water  monomer. 

Figure  6-3  is  a plot  of  the  forces  for  dissociation  of  the  silica  model  used  by  Del 
Bene  et  al.  as  shown  in  Figure  6-1,  where  the  broken  Si-0  bond  is  treated  as  the 
distinguished  coordinate.  The  force  curve  resulting  from  the  transfer  Hamiltonian  is  a 
fairly  good  representation  of  the  second  order  curve.  However,  the  fit  is  not  as  good  as 
seen  for  the  other  systems.  This  is  due  to  the  reduced  number  of  parameters  available  for 
this  parameterization.  Recall  that  only  the  silicon  parameters  are  variable  for  the  fit 
because  those  of  the  oxygen/hydrogen  are  fixed  by  the  water  parameterization. 

Simulations  Of  Hydrolytic  Weakening 
Given  the  transfer  Hamiltonian  parameterization  for  water+silica,  two  molecular 
dynamics  simulations  were  performed  which  demonstrate  the  process  of  hydrolytic 
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Figure  6-3.  Comparison  of  forces  for  dissociation  of  fluorinated  silica  model  using 
CCSD  and  TH-CCSD  potentials. 

weakening.  In  the  first,  a small  SiC>2  ring  was  subjected  to  radial  strain  in  the  presence  of 
the  water  dimer.  The  second  simulation  involved  the  nanorod  which  was  strained 
uniaxially  as  before,  but  with  water  included.  The  stress-strain  curves  resulting  from 
both  simulations  show  a marked  decrease  in  tensile  strength  in  the  presence  of  water 
which  is  an  example  of  hydrolytic  weakening. 

Hydrolytic  Weakening  In  Silica  Ring 

Recent  investigations  indicate  that  small,  highly  strained  two-rings  exist  at  the 
surface  of  amorphous  and  crystalline  silica  [75].  Due  to  high  strain  energy  in  these  rings, 
they  readily  hydrolyze  in  the  presence  of  water  and  are  the  sites  at  which  water  will 
preferentially  attack  the  silica  surface.  Bromley  et  al.  [76]  suggested  a hierarchy  of  silica 
rings  and  chains  as  models  of  extended  silica  surfaces  due  to  the  agreement  between  the 
computed  molecular  ring  vibrational  frequencies  and  IR  bands  measured  on  silica 
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surfaces.  As  an  initial  test  for  the  transfer  Hamiltonian,  a fracture  simulation  of  the  ten 
member  (Si02)io  ring  from  Bromley  et  al.  (Figure  6-4)  in  the  presenece  of  the  water 
dimer  was  performed  using  the  transfer  Hamiltonian.  In  this  simulation,  strain  was 
imposed  by  prescribing  a fixed  radial  velocity  to  each  silicon  atom  and  allowing  the 
remaining  atoms  to  relax  at  each  time  step. 


Figure  6-4.  Ten  member  silica  ring  with  water  dimer.  The  large  circles  are  silicon  and 
the  smaller  are  oxygen  and  hydrogen  atoms. 

Snapshots  of  the  initial,  intermediate,  and  final  configurations  of  the  simulation  are 

shown  in  Figure  6-5  where  it  is  seen  that  the  dual  proton  transfer  mechanism  of  Del  Bene 

et  al.  does  indeed  occur  in  the  simulation.  The  effect  of  the  water  dimer  on  the  strength 

of  the  ring  can  be  observed  in  Figure  6-6  which  is  a plot  of  the  radial  component  of  the 

force  summed  over  the  ten  silicon  atoms.  This  type  of  plot  is  analogous  to  a stress-strain 

curve  which  would  require  division  by  a cross  sectional  area  normal  to  the  loading 

direction  which  does  not  exist  for  the  radially  stressed  silica  ring.  The  strength  of  the 

molecular  ring  is  smaller  in  the  presence  of  water  due  to  the  water  assisted  rupture  of  the 

Si-O  bond  of  the  ring  which  exemplifies  hydrolytic  weakening. 
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Figure  6-5.  Initial,  intermediate,  and  final  stages  of  the  simulation. 

Hydrolytic  Weakening  in  Nanorod 

The  nanorod  introduced  in  Chapter  five  was  also  used  to  demonstrate  hydrolytic 
weakening.  In  this  simulation,  a defect  notch  was  placed  in  the  nanorod  by  removal  of  an 
oxygen  atom  as  shown  in  Figure  6-7.  This  was  done  so  that  the  water  dimer  would  have 
a preferential  site  of  attack  when  in  the  vicinity  of  the  defect  and,  as  before,  uniaxial 
strain  was  applied  to  the  15  atoms  in  the  end  caps  at  the  top  and  bottom  of  the  nanorod. 
This  simulation  contained  113  quantum  mechanical  atoms  with  997  basis  functions  in 
total  and  required  approximately  6 days  of  computer  time  to  complete. 

Snapshots  of  the  initial,  intermediate,  and  final  stages  of  the  simulation  are  shown 
in  Figure  6-8  and,  as  for  the  molecular  ring,  the  dual  proton  transfer  mechanism  is  seen  to 
occur.  A plot  of  the  stress-strain  curves  for  the  full,  notched,  and  notched+water  nanorods 
is  shown  in  Figure  6-9.  There,  it  is  seen  that  the  full  rod  is  weakened  by  the  introduction 
of  the  defect  and  even  moreso  in  presence  of  water. 

Conclusion 

This  chapter  discussed  the  phenomenon  known  as  hydrolytic  weakening  where  the 
strength  of  silica  is  reduced  in  the  presence  of  water.  A proposed  mechanism  for  the 
process  was  outlined  and  a transfer  Hamiltonian  parameterized  which  modeled  the  ab 
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initio  forces  for  the  water  monomer/dimer,  and  silica  in  the  presence  of  water.  Two 
exemplary  simulations  were  presented  which  demonstrate  the  proton  transfer  mechanism 
and  show  that  the  strength  of  the  silica  models  used  is  diminished  when  water  is  present. 
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Figure  6-6.  Plot  of  radial  stress  for  ring  and  ring+water. 
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Figure  6-7.  Nothced  nanorod. 


77 


Figure  6-8.  Initial,  intermediate,  and  final  stages  of  nanorod+water  simulation. 
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Figure  6-9.  Stress  strain  curves  for  nanorod  and  nanorod+water  dimer. 


APPENDIX 

CONVERGED  PARAMETER  VALUES 


Table  A-l.  AMI  core  repulsion  parameters  for  ethane  transfer  Hamiltonian. 


FN1 1 

H 

0.119189531 

FN21 

H 

5.000296781 

FN31 

H 

1.187005328 

FN12 

H 

-5.30830078771 55D-02 

FN22 

H 

5.000124589 

FN32 

H 

1.804518427 

FN13 

H 

-9.0847048920079D-02 

FN23 

H 

2.000275061 

FN33 

H 

2.099731286 

FN1 1 

C 

7.3065280087242D-02 

FN21 

C 

4.997496108 

FN31 

C 

1 .54935964 

FN12 

C 

3.4458410577681  D-02 

FN22 

C 

4.999451186 

FN32 

C 

1.881061799 

FN13 

C 

-3.76078454861 72D-02 

FN23 

C 

4.994770304 

FN33 

C 

2.088195795 

FN14 

C 

1.1 448481 21 801 3D-02 

FN24 

C 

5.001200215 

FN34 

C 

2.682312427 
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Table  A-2.  AMI  core  repulsion  parameters  for  silane  transfer  Hamiltonian. 


FN1 1 

H 

0.137564937 

FN21 

H 

5.000159773 

FN31 

H 

1.196210001 

FN12 

H 

-9.4826688340866D-04 

FN22 

H 

5.000003974 

FN32 

H 

1.798502571 

FN13 

H 

-3.32022091 51 372D-02 

FN23 

H 

2.000166665 

FN33 

H 

2.097608564 

FN1 1 

SI 

0.2578874 

FN21 

SI 

8.99972034 

FN31 

SI 

0.925522328 

FN12 

SI 

6.1 808735024902D-02 

FN22 

SI 

5.009880484 

FN32 

SI 

1.675096047 

FN13 

SI 

-7.20052551 9731 3D-02 

FN23 

SI 

5.000113621 

FN33 

SI 

3.194398939 
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Table  A-3.  AMI  core  repulsion  parameters  for  pyrosilicic  acid/silica  transfer 
Hamiltonian. 


FN1 1 

SI 

0.244008364 

FN21 

SI 

9.001939022 

FN31 

SI 

0.864725617 

FN12 

SI 

4.1 29271 1037864D-02 

FN22 

SI 

5.000018741 

FN32 

SI 

1.949114673 

FN13 

SI 

-4.809251 29691 94D-04 

FN23 

SI 

4.999732104 

FN33 

SI 

2.972359776 

FN1 1 

0 

0.217315862 

FN21 

0 

4.978545859 

FN31 

0 

0.818795183 

FN12 

0 

4.271 61 001 92577D-02 

FN22 

0 

6.999606981 

FN32 

0 

1 .526792547 
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Table  A-4.  MNDO/d  + AMI  core 
transfer  Hamiltonian. 

UDD 

BETAD 

ZD 

ZSN 

ZPN 

ZDN 

UPP 

BETAP 

ZP 

GSP 

GPP 

GP2 

HSP 

ALPHA 

ALPHA 

A1 

B1 

Cl 

A2 

B2 

C2 

A1 

B1 

Cl 

A2 

B2 

C2 


parameters  for  water  monomer/dimer 


-15.67446636 
-10.49396769 
0.932435114 
1 .2567751 1 
1.244652126 
1 .934697891 
-1.567965104 
-0.270500345 
7.88E-02 
7.327374489 
3.626974692 
2.820302283 
1.186342464 
3.278924663 
2.867667176 
-0.786143499 
6.182851647 
0.951527347 
1.255861236 
8.633959981 
2.244723778 
1.032905318 
4.805637752 
1 .960955941 
-0.947660325 
4.419691962 
2.03102124 


repulsion 

8 

8 

8 

8 

8 

8 

1 

1 

1 

1 

1 

1 

1 

8 

1 

8 

8 

8 

8 

8 

8 

1 

1 

1 

1 

1 
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